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I . Introduction 


It is expected that large space structures will be placed into orbit in 
the not-too-distant future. Such structures will lack the damping forces due 
to ground reactions or to hydraulic or aerodynamic forces available to earth- 
bound structures, thus, if they are excited dynamically by docking maneuvers 
or by control reactions, they could be expected to continue vibrating for 
hours or even for days. 

It is this problem which is addressed in this report. The work to be 
discussed divides roughly into two areas. In the first, the questions such 
as how much damping is required, how are the dampers to be located, and how 
their required performance can be specified are examined. In the second, the 
design of the dampers is investigated, and is supported by the demonstration 
of a prototype damper. 

Work on damping requirements resulted in five computer programs and 
proofs of two criteria for the optimum location of dampers. Also, a method 
of minimizing “control spill over” by locating actuators was demonstrated, 
and the characteristic equation of an undamped system using reanalysis theory 
was derived. 

Work on damper design included a review of possible systems. Since 
hydraulic and pneumatic systems were eliminated for poor reliability, only 
electromagnetic or piezo-electric systems were considered. Passive dampers 
were eliminated because they need very fine tuning. Thus, only linear or 
rotational dampers using active electromagnetic or piezoelectric drives were 
left. A prototype linear damper using a moving magnet and a fixed coil was 
built and demonstrated using a current feedback proportional to local 
structure velocity. 



II. Technical Developments 


1. Computer Programs 


Several FORTRAN subroutines and programs were developed under 
this grant. These are documented in Appendix A. The capabilities of 
the programs are enumerated below. 

Subroutines ASMD1 and ASMD2 compute the complex eigenvalues of a 
damped system using the undamped modes of the structure as the 
assumed modes (assumed mode method). In subroutine ASMD1, a diagonal 
damping matrix must be used, whereas in ASMD2, a general damping 
matrix may be used. 

Subroutine SPSTGN solves the damped eigenvalue problem using a 
reanalysis approach. It assumes a diagonal damping matrix. 

Program COPZ computes the optimal (i.e., minimum) gains of a 
diagonal damping matrix for specified damping ratios. 

Program NORMAL is designed to rescale the mode shapes to unit 
generalized mass and to make the rigid bodies orthogonal to each other 
with respect to the mass matrix. 

These programs are stored in the Langley computer. 


2. Proof of Optimal Damper Location 


In previous papers [Ref. 1 and 2] we have proposed two criteria 
for the optimal location of a damper. They are the Minimum Constrained 
Frequency Criterion (MCFC) for a single damper [Ref. 1} and the 
Maximum Frequency Separation Criterion (MFSC) for multiple dampers 
[Ref. 2] : 

MCFC s The optimal damper location is where the constrained 
frequency a^ is a minimum. 

MFSC: The optimal locations of dampers are where the con- 
strained frequency associated with the damper location 
has the largest separation from the corresponding 
undamped natural frequency of the system. 
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In Ref. 1 and 2. these criteria are demonstrated using the 
Langley beam and grillage models. These demonstrations show that in 
choosing between sets of positions of dampers, the one that gives the 
best results is the one that gives the greatest separation between 
the undamped natural frequency and the corresponding constrained 
frequency. In this section, an analytical proof of MCFC for a single 
damper will be presented . The proof of the more general MFSC is 
still under investigation. 

Proof of MCFC 


The proof of MCFC is based on the existence of fixed 
points in the frequency response curve when a single damper is 
introduced into an n-dof undamped system. This is an extension of 
the classical damped vibration absorber theory [Ref. 3], in which 
Den Hartog showed that when a damped vibration absorber is attached 
to a undamped SDOF main system, there are points on the frequency 
response curve which are independent of damping. These points are 
the fixed points [Ref. 3]. In this section, we will first show 
the existence of fixed points when a single damper is introduced 
into an undamped system. The responses at the fixed points can then 
be evaluated from the undamped equation of motion. The MCFC cam 
then be proved using the characteristics of the resonance curve of 
the undamped system near resonance. 

Consider the case of an n-dof undamped system under sinusoidal 
excitation 


[M] {x} + [K] { x } = 


( 1 ) 


where [M], [K] are mass and stiffness matrices and {x} is the 
displacement vector amd {F} is the magnitude vector of the forcing 
function. When a single damper at dof J is introduced into this 
system, the equation of motion becomes 

[M]{x> + [C] {x} + [K]{x} = {F}e jwt (2) 

where 

[C] = c{e J }{e J } T (3) 

{e j* T = [° • • • 0 1 0 ... 0] 

$th component 


Let the steady state solution of ( 2 ) be 


{x} 


{X}e 




(4) 
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Then, use (3) and M) the steady state response {X} can be computed 
from 


[Z(o))]{X} = CF} - C(jj{e J }{e J } T {X} 


( 5 ) 


where 


[Z(w) ] 


- 0 ) [M] + [K] 


Define 


[R(w)l - [Z(a>) ] 
and note that 


( 6 ) 


{e J } T {x} - Xj (7) 

[R] { e j} = (Rj) 

where (R ) is the Jth column of [R], Premultiply ( 5 ) by [R] and use 
J 

( 7 ) to obtain 

{X} = {X q } -jwc Rj Xj (8) 

Note that 

(X q } = [R(w) ] {F } 

= steady state solution of the undamped system. 


The Jth eauation of ( 8 ) is 


X -jwc R X T 
oJ JJ J 


or 


X. = 


QJ 

l + jwc r 


JJ 


(9) 


Equation ( 8 ) together with ( 9 ) are the general frequency response 
reanalysis equations when a damper c is introduced at dof J. It is 
interesting to note that at the frequency a at which 


R JJ ( V ■ 0 


(9a) 
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the response x is independent of damper gain c. Thus, we have 

shown the existence of fixed points for the response curve X (w). 

J 

Furthermore, those frequencies are the antiresonant frequencies for 

2 

dof J, A typical frequency response curve for X (or ai X T> the 

J J 

acceleration response) is shown in Fig. 2.1. In Fig. 2.1, point A is 
the fixed point corresponding to the smallest root of Eq. (9a). By 
definition of fixed point, all response curves pass through A 



to 


Fig. 2.1 Frequency Response Curve for [x^j Showing fixed point A 
a is the frequency corresponding to point A. 



regardless of the damper gain c. Thus, it is obvious the optimal c 
< in terms of the smallest response ) should produce a response curve 
that peaks at point A. In other words, the optimal response curve 
has a slope zero at the fixed point A. 

Prom the above discussion, we know that once a damper location 
is selected, the minimum response X is determined by the undamped 

response at the fixed point frequency. Since near resonance, the 
undamped response is monotonically increasing as the frequency is 
approaching the undamped natural frequency from below, (see Pig. 2.2 ), 
we conclude that the optimal location for lowest amplification for a 
particular mode with undamped natural frequency 0 ) is where (<*) -a ) 
is a maximum. n n 

This proves MCFC for minimum frequency response. Our original 
MCFC pertains to the modal damping ratio. Since modal damping is 
closely related to the amplification factor we conclude that MCFC will 
produce a design with maximum damping ratio. 

Discussion 


(1) In general, there are ( n-1 ) fixed points for response X j( where 
n is the number of elastic modes of the system. 

( 2 ) The responses X^ (I?<J) exhibits fixed points at frequencies 
where 


,X r 1 n " ,X t 1 

I c=0 I c=°° 


or from ( 8 ) and ( 9 ) , condition ( 10 ) becomes 


lx 0.' - IX 0! - R IJ < " )X OJ l 


which corresponds to the following two equations 


( 10 ) 


(ID 


and 


OI 


= X 


01 




01 


-X 


01 


R T .Jw)X rt , 
I J 0J 


(12a) 


(12b) 


(3) Our proof only pertains to X^., however, experience showed that 
once Xj(w) is minimized, the other XjfuO's are almost optimal. 

(4) The proof cannot be extended to multiple dampers since no fixed 
points exist for the multiple damper case. 
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for mode n 


Fig. 2.2 Frequency response of undamped n-dof system. 
jXj (u)| increases monotonically as 
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3. Locating Actuators for Minimum Spill - Over 


A frequently discussed problem in the control of LSS is that of 
"control spill-over" • A method of minimizing the effect of control 
spill-over is proposed in this section. It is shown that to control 
modes using actuators, the spill over effect can be minimized 

by insuring that n secondary modes are not excited. This requires 

s 

n > n . The method described allows only n* = n - n of the n 
as a s a 

actuators to be controlled independently. Conditions of actuator 

placement is also indicated. 

Consider a proportionally damped system whose finite element 
model can be written as 

Mx + Cx + Kx = Bu 


where M, C, and K are mass, damping, and stiffness matrices of the 
system respectively, and x is the system displacement vector, 
u is the control vector and B is the actuator influence matrix. 

In modal coordinates, (13) becomes 
q + 2Ewg + wq=^Bu 


where <t> is the modal matrix and is the solution of the undamped 
eigenvalue problem 

2 

w, M p . =* K p . 
i i i 


E = <5 C <X> 


n 


U) 


w. 


w. 


GJ 


n 


In general, only a few of the important modes are controlled. 
Designate these modes as <J c and the remaining modes as 
q^, the residue modes. Then Eq. (14) can be written as 


(13) 


(14) 


(15) 
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( 16 ) 


• 2 T 

q +2Eo>q + u q =4> Bu 

c c c c c c c 


and 


(17) 


. 2 T 

q + 2E u q +cj q = <£ Bu 

r r r r r r r 


How the physical response is given by 


x = <t> q + <5> q 
cc rr 


(18) 


Hence, if q 0, the physical response is influenced by the 
uncontrolled mode. This effect is called control spill-over . 

T 

Examination of ( 17 ) shows that / o if B u 0. 


It would appear then if we choose B such that 


B - 0 
-r- 


(19) 


the control spill over problem is solved. Unfortunately, condition 
( 19 ) can not be met in general . 

In the following, a simple technique of reducing the control 
spill-over effect is described. Partition the residue modes into 
a set of secondary modes p and higher modes, with the higher modes 
having negligible effects <?n the system. In this way, we have 

x=4>q+<X>q (20) 

-c-c -s-s 

How, to minimize control spill-over, minimize q , or ultimately 
make q g = 0. The response q g can be solved from 

2 T 

q+2Ewq+(jq=<t>Bu (21) 

How if 

<I> u = 0 (22) 

s 

then 



Rewrite (22) as 
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A u =0 where A = <t> B 

- - 8 


( 23 ) 


Assume n < n , then partition A and u into 

& B. “ ” 


A = [A^ • A^] f with A^ nonsingular 


u 


u , 

l 


u 


(24) 


Using (24), (23) becomes 


A.u. + A,u_ = 0 
li d d 


(25) 


From ( 25 ) , solve for 


Thus, if 


u 


- 

-1 


3 -A _ 

A. 

u . 

-d 

-l 

-l 

u. 1 


I 

- 1 




a 




-1 

u . 


-A , A. 

-d 


d l 


u . 

-l 


(26) 


(27) 


Then the condition ( 22 ) will be met . Check 


4> B u = A u = [A. A ] 
s - — - — l a 


-A ~ 1 A . 
d l 


u . 

-l 


= (A . - A . )u . =0 
l ll 


Equation (25) implies that out of the n actuators, only n' 
n of them can be controlled independen?ly . Define 


n - 


D 


I 

-A ^A . 
d l 


(28) 


then 


u - D u . 

- -l 


( 29 ) 
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Substitute ( 29 ) into ( 16 ) to get 

• 2 T 

q+2E(jq+wq“<I> B D u. 
c c cc c c -c - - -l 

n xn nxn n xn' n'xl 
c a a 


The design problem is then to choose u^ to control 


in 


where 


2 T 

q +2 Ecj q + w q = «J> B'u' 
c c c c c c -c 


B' - B D = B 


u = u . 

- -i 


-a” 1 A. 
d l 


A - [O B] 
-a 


Mote that the placement on actuator should be such that one can 

find a non-sinqular n xn matrix A. from the matrix A. 

s s — d — 


4. Characteristic Equations for Undamped System with 
Generalized Control Law 


The characteristic equation of an undamped system with a 
generalized control law is derived using reanalysis theory in this 
section. This equation can be implemented in the computer programs 
for efficient eigenvalue analysis or control gain synthesis. 

Consider an undamped n-dof system with n -actuators, the 
equation of motion can be written as 

Mx+Kx=Bu 


where M. K are the nxn mass and stiffness matrices 
respectively, u is a n & xl vector, i.e.. 



U 

n 


11 1 


(30) 


(31) 


(32) 


u = 


B is the nxn Boolean distribution matrix with n nonzero 
rows. We will designate and X g as the dof locations 
for actuator and sensors respectively. Now define 

T 

X = B x = n xl vector 
-a -a - a 

(33) 

T 

X = B x - n xl vector 
-s -s - s 

then we conclude that 


B - B (34) 

d 

Furthermore, assume X is contained in X , and partition 

a — s 

X as 


x 

-a 


x 

-s 



Now, the system dof can be partitioned either as 






* ' 

X = 

X 

-a 

or 

X * 

X 

-a 

X, 


X 

-r 



-b 

X 

-c 


Note that 


x 

-r 


Also note that 


n = n + n = n + n. + n 
a r a b c 


Now, express 


x - B x + B.x 
-s -a- -b- 


Assume the control law is 



(35) 


( 36 ) 
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( 37 ) 


u = - ( c x + c, x, + k x + k, x.) 
-a-a -b-b -a-a -b-b 


Substitute ( 6 ) into ( 1 ) , so that for the closed loop system 

M x + K x = -B (c x + c. x, + k x + k. x, ) (38) 

- - - - -a -a -a -b -b -a - -b -b 

nxn n xn, n xn, n xn, n xn L 

aab ab ab ab 

or, assuming 


x = Xe 


st 


(39) 


2 -1 

and letting R = (s M + K) equation (38) becomes; 


x => -R B (sc x + sc u x. + k x + k. x. ) 
- -a -a -a -b -b -a -a -b -b 


or 


x > -R B [(ac + k)x + (sc. + k. )x. ) { 

- -a -a -a -a -b -b -b 


(40) 


and finally 


x = -R B [ ( sc + k ) (sc. + k, )] 
~ - -a —a -a ~ b - b 


x 

-a 


^b 


Now, premultiply (40) by B to produce (41) and premultiply 
T ~ a 

( 40 ) by B. to produce ( 42 ) . 

— b 


T T 

B x = -B R B [H H, ] 
-a - -a - -a -a -b 


x 

-a 


L -b 


(41) 


T T 

B x = -B w R B [H H. ] 

-b - -b - -a -a -b 


x 

-a 


^b 


(42) 


Use the definition 


x = B x 
-a -a - 


*b ■ 5b S 
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n.xl n, xn, n.xn' n’xl 
b bob 


where 


R = portions of R associated with x only 

SIS “*d 


R^ a = portion of R associated with the coupling 


between 


x and 
-a 


^b 


Rewrite ( 43 ) as 


or 


X 


-R 

H 

a 


aa 


x. 


-R, 

H 

b 


ba 


A 

x 


A 

X = 


-R 

aa 


H 



A 

X 


det[I + 


R H 

aa - 



0 


thus 


det 


+ R H 

aa -a 


R 

aa 



R. H 
ba -a 


[I + 


R , H. ] 
ba -b 


0 


( 44 ) 


14 1 


d 



is the characteristic equation of the closed loop system where 


H “SC + k 
-a -a -a 


-b 


3 ^b + Sb 


Example 

Por the 2-dof system shown in Figure (4.1) 


°1 “ ~ (C a X l + C b X 2> 


H = sc 
a a 


H. = SC, X = X, 

b b -a 1 


^b = X 2 


R = R-. 

aa 11 


R bb = R 22 


R ~ R_ , 
ba 21 


Thus, Eq. (14) becomes 


det 


1 + 3c a R 11 
3c a R 21 


3C b R ll 

3C b R 21 


- 0 


-> X, 


U i > 



Fig. 4.1 
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5 . Damper Des ign 


Possible methods of controlling vibrations on large space 
structures have been reviewed. Two basic approaches are identified: 

1. Modification of the control system so that the poles are well 
damped. This approach is based on the assumption that the 
major, if not only, cause of vibrations is the control 
system itself. 

2. Addition of a separate damping system, which might be: 

2.1 Passive damping 

2 . 2 Active damping 

The working medium of a ground based damper of either kind 
might be electromagnetic, but it would more probably be hydraulic 
or pneumatic, or, in the case of passive systems, magnetic or 
frictional. For space applications, hydraulic, pneumatic, and 
frictional systems must be ruled out as too unreliable in the space 
environment when little or no maintenance is possible. This 
leaves electromagnetic, magnetic, and possibly new systems based 
on the piezoelectric or other effects. 

A list of possible systems is given in Table 5.1. The active 
damping systems using electromagnetics are assumed to incorporate 
coils moving in the fields of permanent magnets. The current in a 
coil would be provided by an electrical circuit in response to the 
output from one or more sensors. For the purpose of this study, such 
a circuit is assumed to introduce a negligible weight penalty, whereas 
the permanent magnet is considered to be a significant weight item. 
Similarly, the permanent magnet in the eddy current damper of a 
passive system is also considered to be a significant weight item. 

For systems which require a moving mass, such as the seismic mass of a 
passive system, or the mass in an active inertial system, it would 
therefore seem to be advantageous to combine the roles of moving mass 
and magnet in order to save weight. On this basis, the gyro concept 
might have a considerable weight disadvantage, because it would be 
impossible to combine the magnet and the flywheel. 

Table 5.1. Possible Damping Systems 


1.0 Modification of Control System 

2 . 0 Damping System 

2.1 Passive Damping: Magnet and conductive strip 


2.1.1 

Seismic mass with stationary magnet 

2.1.2 

Magnet as seismic mass 

2.1.3 

Gyro 

2.1.4 

Two-force member 

2.2 Active 

Damping: Electromagnetic 

2.2.1 

Separate mass and magnet 

2 . 2 . 1 . 1 Linear 

2. 2. 1.2 Rotating (inertia wheel) 

2.2.2 

Moving magnet 

2 . 2 . 2 . 1 Linkage system 

2 . 2 . 2 . 2 Rack and pinion 

2.2.3 

Gyro 

2.2.4 

Two-force member 

2 . 3 Active 

Damping: Piezoelectric 
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In comparing active against passive damping, the relative 
disadvantage of requiring an electrical system for the former must 
be weighed against the need to tune the seismic mass to a given 
frequency, so that the latter is essentially a narrow band device, 
whereas the active damper is relatively broadband. 

The systems mentioned in Table 5 . l are discussed in more detail 
below. 

Modification of Control System (1.0) 

Modern control theory prescribes methods for designing 
control circuits having poles in any desired location, thus active 
damping can be achieved directly at the major source of disturbances. 
A very simple example is illustrated in Figure (5.1) which shows a 
one-dimensional system consisting of two masses connected by a 
spring, having one rigid-body translational mode, and one vibrational 
mode. Suppose that the impulse from a thruster is applied at one of 
the two masses . The vibrational mode will be excited and the system 
will move off with one mass coming to rest every half cycle. Now 
suppose that there are two thruster impulses applied one half-period 
apart, then the system will move without any internal vibration after 
the first half-period. Applied on a very much larger and more 
complex scale, such an approach could be used to maneuver a large 
space structure without any residual vibrations of serious amplitude. 

Passive Damping (2.1) 

The only reliable method of passive damping appears to be 
through the use of a magnet and a conductive strip. When the strip 
is passed between the poles of a magnet, eddy currents are induced in 
a direction normal to the velocity and to the magnetic field, and the 
effectiveness of such a damper is directly proportional to the 
product of electrical conductivity, magnetic field strength, and pole 
area. Unfortunately, the damping forces so produced are relatively 
small at frequencies of interest. 

Four possible systems are shown in Figures (5.2) to (5.5). The 
first, system 2.1.1 in Figure (5.2), is a seismic damper, while for 
system 2.1.2 in Figure (5.3) the magnet doubles as the seismic mass . 
For system 2.1.3 in Figure (5.4), the magnetic damper is used in 
conjunction with a gyro to provide damping against angular motion of 
the structure. Finally, for system 2.1.4 in Figure (5.5), the 
magnetic damper is used in a two-force member spanning part of the 
large structure. 

Active Damping: Electromagnetic (2.2) 

In the active damping system, energy is removed from the 
structure by the action of an electromagnetic actuator or an electric 
motor. For the purpose of illustration, a voice-coil shaker of the 
type used for vibration testing has been assumed for the analyses of 
the active systems in the appendices . Such shakers are relatively 
heavy, because they contain powerful permanent magnets, although 
careful redesign for space applications might result in considerable 
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DOUBLE IMPULSE, 2nd IMPULSE APPLIED 
ONE HALF PERIOD AFTER 1st IMPULSE 

Fig. (5.1) 1.0 Modification of Control System 
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weight reductions. However, the weight penalty is directly 
proportional to the maximum energy Which can be removed per cycle, as 
is shown in Appendix B. The design problem with a damper of this 
kind is that of ensuring that this maximum energy reduction per cycle 
is in fact available. 

The system would be driven through a current amplifier from 
the difference between two integrated accelerometer signals, one 
from the structure, and one from the moving mass. Thus a force 
would be applied to the structure which would be directly 
proportional to the local velocity of the structure. 

The major problem to be faced is how to provide an appropriate 
reaction on the structure so that as much energy is removed per cycle 
as is possible. One approach, shown as system 2. 2.1.1 in Figure 
(5.6), employs a moving mass. As is shown in Appendix B, this system 
can be designed for peak performance at a given frequency. Above this 
frequency, performance is limited by the maximum available shaker 
force, while below this frequency, performance is limited by the 
maximum displacement available. An analogous rotarysystera, 

2. 2. 1.2, shown in Figure (5.7) uses an electric motor driving an 
inertia wheel. 

System 2. 2. 2.1 in Figure (5.8) is an improvement on the previous 
system, in which the shaker magnet doubles as the moving mas3. A 
second version is shown in Figure (5.9). In Figure (5.8), the 
actuator motion is amplified by a "lazy tong" linkage. Such a linkage 
would not be practical, but could be used in a laboratory 
demonstration. In the second system. 2. 2. 2. 2., Figure (5.9). the 
shaker is replaced by an electric motor which drives itself up and 
down a shaft with a rack-and-pinion gear. 

System 2.2.3 shown in Figure (5.10) is similar to system 2.1.3 
3hown in Figure (5.4) except that the passive damper has been replaced 
by an active electromagnetic damper , while system 2.2.4 shown in 
Figure (5.11) is a two force member corresponding in a similar way 
with system 2.1.4 in Figure (5.5). 

Active Damping: Piezoelectric (2.3) 

Acoustical transducers have been made from polyvinylidene fluoride 
(PVF z ) sheet, aluminized on both sides. The sheet is polarized, so 
that when a voltage differential is applied between the aluminized 
coatings, the material strains in one direction. Conversely, when it 
is strained, a voltage differential is induced. If the material is 
bonded to the surface of a structural element, a surface shearing 
force can be induced which will load the structure , and , if properly 
controlled, is capable of introducing damping. Proper use of such a 
material evidently depends on new approaches to structural dynamic 
analysis, but is certainly worthy of consideration. 

Electronics 

The active systems would be driven by feedback circuits. Such 
a circuit is shown in Figure (5.12), which repeats system 2. 2^1. 2 ^f 
Figure (5.8) with the addition of a circuit which compares md y/dt and 
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Fig. (5.11) System 2.2.4 Two Force Member Active Damper 
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FIG. 12 SYSTEM 2.2. 2.1 ACTIVE DAMPER Showing possible control loop 

added to system to yield desired characteristics 










cdx/dt as derived from two accelerometers, the difference being used 
to drive the voice-coil current through a amplifier. 

Comparative Performance 

For a preliminary comparison of performance capability, two 
figures of merit can be considered. 

1. Energy removed per cycle per unit mass at design 

freouency, W /in 
D T 

2 . Bandwidth 

In Appendix B, the value for W D /® T achieved with system 2. 2. 2.1 in 

Figure (5.8) under design conditions is shown to be .0619 J/kg. which 
is 79% of the maximum possible capability of a voice coil shaker. 

This was twice as good as the best system 2.2.1 arrangement of Figure 
(5.6). Therefore, it is tentatively concluded that other active 
damping systems, such as the inertia wheel, in Fig. (5.7). or the two 
force member, in Fig. (5.11), cannot be much better if they are based 
on electromagnetic action. Passive systems, being based on magnets, 
are probably no better, although this has not been investigated. 

Performance of the active system, based on rate of energy 
removal, is constant above the design frequency, but rolls off below 
it. In contrast, the passive system is narrow-band and has to be 
tuned to a given frequency. 

It is shown in Appendix C that the FVF piezo-electric 

material has the potential for work per cycle of as much as 4 J/kg, 

which is more than 30 times that of the voice-coil shaker. On the 

other hand, using the minimum property values available, the work 

per cycle could be as low as .005 J/kg, which is less than one 

tenth as good as a voice-coil shaker. One of the most critical 

properties is the dielectric field strength, E . since a dielectric 

w 

breakdown causes an explosive failure, it might be necessary to 
allow quite a large safety margin, so that PVF 2 might prove no 
better than voice-coil shakers. 

Conclusion of Review 

It is concluded from this series that an electromagnetic 
actuator, suspended in such a way that the magnet forms part of the 
moving mass, is the most promising approach for immediate 
application . 

Design of Experimental Damper 

It was decided, as a result of the review, that an electro- 
magnetic damper should be investigated experimentally, and that 
one should be developed for use on the scale model grillage at 
the NASA Langley Laboratory. The following design criteria were 
used initially 

Maximum Force = 1 lbf = 4.45 N 
Design Amplitude = 0.5 in. = 0.0127 m. 

Design Frequency = 0.36 Hz 
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Instrumentation 


Accelerometer on moving mass - NASA supplied 
Accelerometer on grillage - NASA supplied 
Electronics - off grillage - NASA supplied 

A design which meets these requirements is shown in Fig. (5.13 ). 
Description of Damper 

The magnetic flux is supplied by the two toroidal samarium- 
cobalt magnets. The remainder of the magnetic circuit is temporarily 
completed with mild steel, which has a high saturation flux, and can 
therefore be designed for minimum weight. However, pure iron will 
be substituted when it becomes available. The annular gap is large 
to accomodate the windings of the fixed coil, this dictates a design 
which minimizes flux leakage at the gap. For example, two magnets 
are used where one would suffice if the gap were small, also, the 
internal diameter is larger than is required merely for clearance. 

The entire magnet assembly moves along a hardened steel shaft 
on linear bearings, with one inch of useful travel. An 
accelerometer is attached to the moving magnet to feed back the 
acceleration of its mass, which is directly proportional to the 
damping force. The fixed coil has a Delrin core, and is designed to 
take ten layers of 26 gage magnet wire. 

An outer case of polycarbonate tubing is used to keep the 
linear bearings clean. It also holds a track in which a small ball 
race moves to prevent the magnet from turning and thereby twisting 
the leads to the accelerometer. The end cap of the cover is 
removable so that the accelerometer leads can be connected to 
terminals . 

Magnet Analysis 


The samarium cobalt magnets have the following specifications : 

Outside Diameter = 0.75 in. 

Inside Diameter = 0.43 in. 

Magnetic Induction = 8000 Gauss 

=0.8 Tesla 


The following is calculated: 

Magnetic Flux = 153 micromaxwells 

The remaining circuit is designed so that the magnetic flux 
density in the mild steel does not exceed 66% of the 2.1 Tesla 
saturation flux density. 

For the purpose of the subsequent analyses, the total 153 uM 
flux is assumed to cut the coil windings. 
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Fig. (5.13) N,A,S.A, Active Damper Actuator Prototype 







Analvsis of Coil 


It is assumed that each of the ten layers of 26 gage wire 
contains 67 turns, based on the following analysis 

Spool length = 1.25 in. 

Diameter of bare wire = 0.015945 ins 
Diamater of enamelled wire (measured) = 0.01772 ins 
Winding efficiency =* 67 x 0.01772/1.25 = 0.95 = 95% 
Turns per meter = 10 x 67 x 39.37/1.25 = 21100 m -1 

The relationship between the force acting, F, the flux <J>, the 
turns per meter, n, and the coil current A is developed as follows. 
Assume an effective coil diameter D, and a gap width w , then the 
flux in the gap is g 


B =■ 4>/t t D w 

(5.1) 

gap g 

The length of wire immersed in the gap is 


(5.2) 


L = n D n w 
gap g 


and the force acting between the coil and the magnet is 


(5.3) 


F = I B L - 4>nl 

gap gap 

For a force of one lbf, or 4.45 N 


I = F/n 4> = (4.45)/(21100) (153 x 10 6 ) = 1.38 A 


In operation, this would be a peak value, the actual DC current 
recommended for 26 gage wire is 0.51 A, based on 254 circular mils 
(in) at 500 per A. Thus a damper operating continuously at maximum 
amplitude would be limited to a maximum force of 


F = V? (0.51/1.38) = 0.52 lb 
max 

However, during tests, peak amplitude will vary, so that the 
damping constant can be set to correspond to a much higher maximum 
force without causing excessive heating in the coil. The resistance 
of the coil, based on 0.0410 ohms/ft., and a mean diameter of 1.2 
ins , is 


R = (0.0410) (10) (67) (1.2) v / ( 12 ) - 8.63 Ohms 


Other estimates, based on different values for resistivity, 
have been as high as 23 Ohms. Based on 8.63 Ohms, the peak voltage is 
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Vmax = (1.38) (8.63) 


11.9 V 


Vmin = -11.9 V 


Analysis of Damping Constanta 


If we redefine K so that it is the ratio of structural 
amplitude Ix q I to damper amplitude, i.e.. 


K - IX I/D 
o m 

then, from Eq. (B.2) of Appendix B 

. /v 2 ,^ 2.22 

1/K = 1 + c /u 


also, with the linkage ratio R set to unity, Eq. (B.l) can be 
written as 


(5.4) 


(5.5) 


F/F - wc D K/F (5.6) 

m mm 

For the present design, the following values are anticipated 


m 2 = 0.428 kg 
D - 1/2 inch = 0.0127 m 

ID 

F =1 lbf =» 4.45 N 
m 


From Eq . (5.5), we can construct the following table for the damping 
constants , c . 




Table 5.2 

Values of Damping Constants 


K 

'V 

C/(J 

f = . 35Hz 

C (Ns/m) 
f = 1 Hz 

f = 3.5Hz 

.8 

(m) 

.0102 

(kg) 

.321 

.7059 

2.017 

7.059 

.5 

.0064 

.7413 

1.630 

4.658 

16.3 

.3 

.0038 

1.361 

2.993 

8.551 

29.93 

.1 

.0013 

4.259 

9.37 

26.76 

93.7 


and from Eq. (5.6) we can construct the following table for the 

corresponding force ratios, F/F . 

N 


Table 5.3 Force Ratios 


K 

,X o' 

<m) 

f = . 35Hz 

F/F 

f = 1 Hz 

f = 3.5Hz 

.8 

.0102 

.00354 

.0289 

.354 

.5 

.0064 

.00511 

.0418 

.511 

.3 

.0038 

.00563 

.046 

.563 

.1 

.0013 

.00588 

.0480 

.588 


Note that the full force capability of the actuator is not 
required until the frequency exceeds 3.5 Hz. 

In the proposed design, the damping constant c will be set by 
the electronic circuit. Suppose we pick a value of 2 Ns/m and 
examine the corresponding values for K and F/F^ as given by Eqs . ( 5 . 5 > 
and (5.6). Note that a value for K of unity could correspond to a 
peak-to-peak structural amplitude of one inch while a value for F/F 
of unity would correspond to a maximum force of one lbf, the design 
value. The results are given in Table 5.4. 

Table 5.4 Values for K and F/F when C = 2 Ns/m, (Eqns. 2 and 3) 

n 


f( Hz) 

K 

F/F„ 

.1 

.133 

.00048 

.2 

.261 

.00185 

.3 

.374 

.00403 

.35 

.426 

.00535 

.4 

.474 

.00678 

.5 

.558 

.01006 

.8 

.732 

.0210 

1.0 

.802 

.0288 

1.5 

.899 

0.0484 

2.0 

.937 

.0672 

2.5 

.958 

.0859 

3.0 

.971 

.1045 

5.0 

.989 

.1773 


Note that the maximum force is not approached even at 5 Hz. 

Had a larger damping constant been used, these forces would have 
been larger, but the value for K would have been smaller. For 
example. Table 5.5 shows how K varies with c at frequencies of 0.35, 
1 . 0 and 3.5 Hz . 
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Table 5.5 K vs. c (Eqn. 2) 


c 


Ns/m 

0.35 HZ 

0.1 

.994 

0.2 

.978 

0.3 

.953 

0.4 

.920 

0.6 

.843 

0.8 

.762 

1.0 

.685 

1.5 

.532 

2.0 

.426 

2.5 

.352 

3.0 

.299 

5.0 

.185 

10.0 

.094 


K 


0 HZ 

3.5 Hz 

999 

1.000 

997 

1.000 

994 

1.000 

989 

.999 

976 

.998 

958 

.996 

937 

.994 

873 

.988 

802 

.978 

732 

.966 

607 

.953 

474 

.883 

260 

.685 


None of the examples shown in Table 5 . 5 correspond to a value 

for P/F exceeding unity. Thus, the present design, if achieved, 
m 

should be adequate for all of the cases shown, provided that c 
is set at 2 Ns/m or less. 

Preliminary Tests 

A prototype unit was built from the design shown in Figure 5.13. 
a photograph of this unit is shown in Figure 5.14, while the analog 
circuit used in an initial test of the unit attached to the NASA 
grillage is shown in Figure 5.15. 

Because of the possibility of hitting the stops when the moving 
mass is not centralized, a spring was added to the system which gave 
it a natural frequency of about 2 Hz. Using switch no. 2, the system 
was excited at a natural frequency close to 4 Hz., then it was 
switched to feedback control, and the resulting damping was observed. 
Through switch no. 1, the inner feedback loop incorporating the 
accelerometer mounted on the moving mass could be included at will. 
When it was excluded, rapid damping of the 4 Hz. mode was observed, 
but when the inner loop was included, the system at first damped, 
and then went unstable at a frequency lower than 2 Hz. 

It was concluded that the accelerometer on the inner loop could 
be discarded, but that an inner loop should be incorporated which 
included an LVDT to measure position. This, way, a very low natural 
frequency could be simulated for the system, while a bias voltage 
could be applied to offset gravitational effects. At the same time, 
studies should continue on a system which incorporates an inner loop 
accelerometer . 

Conclusions 

An experimental damper has been designed in which a mass of 
.428 kg ( .95 lbm) moves over a peak-to-peak amplitude of one inch 
under a programmed force which can be as much as one lbf . The 
design includes an accelerometer attached to the mass. 
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Fig. 5.14 Prototype Damper 




SIGNAL 
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Control Circuit for Prototype Damper 















The programmed force is to be produced from the generation of a 
current in an amplifier driven by the difference between a signal 
proportional to the structural velocity and one proportional to the 
acceleration of the mass. Thus a damping force c will be generated 
on the structure. 

Tentatively, the following design values are suggested 

Damping constant = 2 Newton sec/meter 

Design frequency = 0.35 Hz 

Structural Amplitude = .43 inches peak-to-peak 


From early experience with the prototype, it is concluded that 
the inner-loop accelerometer should be replaced by an LVDT if an early 
working system is required. However, research should continue into 
the use of the inner accelerometer, and an alternative to the LVDT 
should be found if dampers with longer strokes are found to be 
desirable . 


III. CONCLUSIONS 


In conclusion, this work has covered two aspects of the problems involved 
in damping large space structures. On the one hand, the analytical problem of 
locating dampers has been investigated, while on the other, the problems of 
damper design have been reviewed. These considerations are summarized below 

(1) Five computer programs have been developed. They are 

(i) asmdi - Assumed mode method with diagonal damping matrix 

( ii ) ASMD2 - Assumed mode method with full damping matrix 

(iii) SPSTGN - Damped eigenvalues using reanalysis 

(iv) COPZ - Optimization of damper gain 

(v) NORMAL - Orthonormatization of mode shapes 

( 2 ) Proof has been presented for two criteria for the optimal location of 
dampers . They are s 

( i ) mcfc - That the optimal single damper location is where the 

constrained frequency is a minimum 

(ii) MFSC - That the optimal locations of several dampers are where 

the constrained frequency associated with the damper 
location has the largest separation from the corresponding 
undamped natural frequency of the system. 

( 3 ) A method of minimizing the effect of control "spill-over" has been 
proposed and demonstrated. 

(4) The characteristic equation of an undamped system with a generalized 
control law has been derived using reanalysis theory. 

(5) A review of possible damper designs has been conducted. It was 
concluded that the most promising designs are active systems, using 
electromagnetic or piezoelectric actuators with linear or rotational motion. 

( 6 ) A prototype linear electromagnetic damper was built and demonstrated 
using a moving permanent magnet. The damper is driven by feedback from an 
accelerometer, mounted in the structure, and integrated to provide a feedback 
force proportional to structural velocity, it was found necessary to 
incorporate a centering spring which will be replaced in future designs by a 
position feedback from an LVDT or other device. 
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Appendix A-l 


Subroutine ASMD1 


Assumed Mode Method - Diagonal Damping Matrix 


I. INTRODUCTION 

asmdi is a FORTRAN subroutine, based on the assumed mode method, 
that can be used to compute the eigenvalues of a damped system. The 
eigenvectors of the undamped system are used as the assumed modes . 
When all modes of the undamped system are used, the method is 
equivalent to a direct solution. Experience shows that retaining 1/3 
to 1/2 of the modes of the original system in the assumed mode method 
usually leads to accurate eigenvalues while providing considerable 
savings in computer time . 

II. ASSUMED MODE FORMULATION 


Let the free vibration of a damped system be described by 

[m]{x> + [c] {x} + [k] {x} = {0} (Al.l) 

where [c] is a diagonal damping matrix. Let {p.}, u>. he the solution 
of the corresponding undamped problem 1 1 

w 2 [m]{p} = [k] {p} (A1.2) 

Furthermore, assume the mode shapes are normalized to unit 
generalized mass, i.e.. 


fp.} T [m]{p.} = 1.0 (A1.3) 

Define 


L 


{X} = 

E {p } q = [4>] 

(q) 

(A1.4) 


* = 1 i i 


nxl 

nxL 

Lxl 



where L is the number of modes used. Substitute (A1.4) into (Al.l) 

T 

and premultiply [$] to get 

P • 

[ I Hq) + [5]{g} + [ A ]{g} = {0} (A1.5) 
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where 


[ I ] = LxL identity matrix 
2 

[ A ] = [ Uj ] = LXL diagonal matrix 
T 

[c] = [<X>] = LxL full matrix 

Set 

(A1.6) 

2Lxl 

Then, (A1.5) becomes 

{ 2 ^} = ~[c][z 1 } - [ A J{z 2 ) (A1.7a) 

By definition (A1.6) 

{z 2 > = { Z]L } (Al.7b) 

Place (A1.7a), (A1.7b) together 



{z} = [A] { z } (A1.8) 

where 


r i r >ri j ■} 



-[c] 

-[ A ] 


-[•] [c][$] 

-t <V i 

[A] - 



= 




[I] 

[0] 


[I] 

(0) 


A standard eigenvalue problem can be formulated from ( 8 ) and solved 
using the routine RESV in the NASA Langley library. 


III. INPUT PARAMETERS 


III.l Formal Parameters 

Description of the 6 formal parameters 

Ns number of DOF's of the original system 
NC: number of dampers 
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IPRTs Dynamical matrix ([A] in (8>) printing flag, 

= O no printing 
= 1 print 

III . 2 Common Blocks 

In addition to the formal parameters, the following variables 
enter the subroutine ASMD through the common block set contained in 
file COMASMD. This file must be placed in the main program. 

NTITLE( I ) s a holerith array of up to 80 columns which is used to 
define the title of the job 
JC( I ) : damper location (DOF) of the Ith damper 

C( I ) : initial damper values 

WN( I ) : undamped natural frequency of the ith mode 

PHI(J,I): mode shape coefficient at dof J of the Ith mode 

NM: number of modes used in the assumed mode method 


IV. OTHER SUBROUTINES USED 

ASMD calls NASA Library routine KESV to compute the eigenvalues 
of the dynamical matrix [A] . The computed complex eigenvalues are 
stored in array ER( I ) and EI(I), which are the real and imaginary 
part of the Ith eigenvalue of matrix [A] . 
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V. USING SUBROUTINE ASMD 


The following is a sample program that calls routine ASMD1 . 


C 

C 


C 

C 


PROGRAM TESTAC1 < INPUT , OUTPUT f 1 APE12 ) 


DIMENSION COO) 

COMMON / DATA1 /C<50) 

COMMON / DATA2 / PHI <?9f99) fUN(??> 

COMMON / DATA3 / EIGDATA ( 40 f 20 f 4 > f DAMP C 40 > 

COMMON / DATA4 / NT1TLEOO) » JC<50> 

COMMON / DAT AS / NM 

COMMON / DATA6 / NDP 

COMMON /DATA7/ER<7?) fEI (7?) 


DO 1 1 3 1 f SO 
1 C(I>=0.0 

C 

C INPUT SECTION 

C 

READ 10»(NTITLL<I)fI~l»00) 
PRINT 10 » ( NTITLE ( I ) f 1 = 1 f 00 ) 
10 FORMAT (00A1 ) 

READ *fNfNMfNC 
PRINT t'N'NM>NC 
DO 20 1=1 f NC 
READ *fJC<I)»C(I) 

PRINT *fJC(I)fC(I> 

20 CONTINUE 

READ * t IPRT 
PRINT * » IPRT 
C 

DO S010 I =1 » NM 
READ* » UNCDfCd) 

PRINT* f UN ( I > f G < I ) 

READ* f (PHKJfI) fJ=1fN) 
PRINT*f(PHI<JfI)fJ*1fN> 
SOU FORMAT (5E1A.7) 

5010 CONTINUE 

PRINT*f *N UN 0* 

DO 9? 1=1fN 

99 PRINT*fI»UN(I)fG(I) 

PRINT*f , PHI* 

DO 999 J=1fN 

999 PRINT*f <PHI< JfI) fI=*1»N) 

C 

C 

C 

C 

CALL ASMDKNfNCfIPRT) 

C 

m - 
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VI. PROGRAM LISTING 


u 

SUBROUTINE ASMD1 ( N » NC t IPRT ) 

C 
C 

C THIS PROGRAM USES THE ASSUMED MODE METHOD TO COMPUTE 

C COMPLEX EIGENVALUES OF A DAMPED SYSTEM. 

C DIAGONAL DAMPING MATRIX. 

C THE PROGRAM IS A MODIFICATION OF PROGRAM CLIO 1/13/C2. 

C 
C 

COMPLEX CMCQF.IG 
COMMON / DAT A1 /C(50> 

COMMON / DAT A2 / PHI < 99 » 9? ) » UN ( 9? ) 

COMMON / DATA3 / EICDA I A ( 40 » 20 » 4 ) » DAMP ( >10 ) 

COMMON / HATA4 / NT I TL E ( 80 ) » JC ( 50 ) 

COMMON / DATA5 / NM 
COMMON / DAT A6 / NDP 
COMMON /DATA7/ER(?9) » E I (99) 

C 

DIMENSION DYEXAC ( 99 » 99 ) » V ( 99 ? 5 ) t UK ( 270 ) 


MAX r 99 
NMD=2*NM 
N = NMD 

FORM DYNAMICAL MATRIX 

NN=2*NM 
DO 100 I - 1 » NN 
DO 100 J = 1 » NN 

00 DYEXAC ( I » J ) -0.0 

DO 200 I - 1 » NM 
I N = I + N M 

DYEXAC! I ? IN) =- UN! I )*UN! I ) 

DYEXAC! IN? I ) ^1 .0 
DO 200 J = 1 ? NM 
DO 200 L- 1 ? NC 
K= JC ! L ) 

DYEXAC! I ? J) -DYEXAC! I ? J > -C < L ) *PHJ (K»I)tPHl (K» J) 

00 CONTINUE 

PRINT*? •*********#****««*»»? OUTPUT 

PRINT 11 

PRINT 10f (NTITLE(I) f I’ It 80) 

0 FORMAT ! 80A 1 ) 

1 FORMAT!//) 

IF! IPRT.EQ.O) GO TO 210 

PRINT 201 

01 FORMAT ! 1H1 »///» # DYNAMICAL MATRIX IS *///) 

DO 202 I-lrNN 

02 PRINT 205 i < DYEXAC < 1 ? J )? J= 1 ? NN ) 

05 FORMAT (/?5X»5E15.6) 

10 CALL RESVIMAXfNfDYEXACfERfEIfOfOfVtUKf 1F.RR) 


•: t * * * * 
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c 


PRINT 1020 » NM 


C 

666 

1020 

N. 

c 

333 

9999 


PRINT *t ****** DAMPERS ****** 

PRINT* » * PLACEMENT(DOF) VALUE* 

DO 666 1 = 1 > NC 

PRINT* » * *»jca)t* *»CQ> 

FORMAT</» 10X,* ASSUMED MODE METHOD » USE *»IS 
PRINT 11 

PR INT* » * THE DAMPED EIGENVALUES:* 

DO 333 I 1=1 » NMD 

PRINT*»II»ER( 1 1 > » “ + * » E I (II)» *J* 

PRINT 11 
PRINT 11 
CONTINUE 

END 


> 


>>* MODES* / ) 
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VII NUMERICAL EXAMPLE 


To illustrate the use of ASMD1, the 2 dof system of Fig. 

( A1 . 1 ) is used . A damper of c = 2,3 is attached tcjSdof 2 . The 
input data and control cards for this sample problem are shown in 
Table 1 . The output is given in Table 2 . 
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Table 1 Control Cards and Input Data for Sample Problem 

The test program TESTAS1 can be run through the use of the 
following job files 


/JOB 

BPU » CM300000 » T 1000 • 

USER » 677121E . 

CHARGE* 102420, LRC. 

DELIVER. 12735 D UVA 
ATTACH, FTNMLIB/UN«LIBRARY. 
GET * TE3TAS1 « 

MAP(QFF) . 

FTN»I 3 TESTAS1 , L 3 0,OPT=1 , B~F. 
LDSET,LIB 3 FTNMLIB. 

F. 

REWIND * OUTPUT • . 

COPYEI, OUTPUT, 0TTAS1. 

REPLACE *0TT AS 1 . 
DAYFILE,L 3 DATAS1. 

REPLACE »D AT AS 1 . 

EXIT. 

REUINDfOUTPUT. 

COPYEI, OUTPUT. X0TTAS1. 
REPLACE. X0TTA31. 
DAYFILE,L 3 XDATAS1. 
PACK.XDATAS1 . 

REPLACE.XDATAS1 . 

/EOR 

test: 2-Dor free system, test 

2 » 2. 1-4 

2 , 2 . 3 * 

0 

0*1 <4 

0,577. 0.577 ^ 

1 . 2247,1 -4 

0.316,-0.400 < 

E0I ENCOUNTERED. 


12735 * UVA 


ASMD 1 - < z 0*1 % ioih M' 0 ; 1 OMfcfi. 
Ofi**frR: oof 2 cr 3 

DfH. f 1 * T. rioT pfirfUQ 


/.llll _/ 


/ 
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Table 2 Output of Sample Problem 


TEST: 2-DOF free system? test ASMD1 


ASSUMED MODE METHOD? USE 


***** DAMPERS ***** 
PLACEMENT (DOF) VALUE 

2 2.3 


THE DAMPED EIGENVALUES: 

1 O.+O.J 

2 -.8999940753023 + 0. J 

3 1 243049 123488+1 . 122804973927 J 

4 -. 12430491 234 88+- 1 . 12280497392/ J 


MUDE 
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Appendix A-2 


Subroutine ASMD2 


Assumed Mode Method - Full Damping Matrix 


I. INTRODUCTION 

ASMD2 is a FORTRAN subroutine, based on the assumed mode method, 
that can be used to compute the eigenvalues of a damped system. The 
eigenvectors of the undamped system are used as the assumed modes . 

When all modes of the undamped system are used, the method is 
equivalent to a direct solution. Experience shows that retaining 1/3 
to 1/2 of the modes of the original system in the assumed mode method 
usually leads to accurate eigenvalues while providing considerable 
savings in computer time. ASMD2 allows a full damping matrix to be 
used. 

II. ASSUMED MODE FORMULATION 

Let the free vibration of a damped system be described by 

[m]{x} + [c] {x} + (kHx) = {0} (A2.1) 

where [c] is a diagonal damping matrix. Let {p.}, w. be the solution 
of the corresponding undamped problem 1 1 

w 2 [m ]{P} =* [k]{p} (A2.2) 

Furthermore, assume the mode shapes are normalized to unit 
generalized mass, i.e., 

{p.} T [tn] {p.} = 1.0 (A2.3) 

Define 


L 


{X} = 

E {p }q = [<X>] 

(g) 

(A2.4) 


* = 1 £ £ 


nxl 

nxL 

Lxl 



where L is the number <^f modes used. Substitute (A2.4) into (A2.1) 
and premultiply by [<t>] to get 

t I Hq) + ic){q) + [ A ]{g) = {0} (A2.5) 

where 
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Set 


[ I ] = LxL identity matrix 
2 

[ A ] » [ Wj ] = LxL diagonal matrix 
T 

[c] = [4>] [c][4>] = LxL full matrix 



2Lxl 


Then , ( A2 . 5 ) becomes 

{z 1 > = -[cHz^ - [ A ]{z 2 ) 
By definition ( A2.6 ) 


(A2.6) 


(A2.7a) 


{z 2 } = { Zl } 

Place (7a), (7b) together 
{z} = [A] {z} 


(A2.7b'> 


(A2.8) 


where 


[A] 


-[c] 

-t a ] 


-[*] tc]fd.] 

-i v J 

[I] 

[0] 


[I] 

[0] 


A standard eigenvalue problem can be formulated from ( A2 . 8 ) and solved 
using the routine RESV in the NASA Langley library. 


III. FORMAL PARAMETERS AND CALLING ASMD2 


The subroutine ASMD2 is called with the statement: 


CALL ASHD2 (NiNM> NROUC » C » PH I » WN » FIR » F.l i NT 1 TLE > IPRT ) 


! 


The formal parameters are defined as follows: 


N= N.D.O.F. IN THE. SYSTEM 

NM= NO. MODES USED IN THE ASSUMED MODE CALC. 

NROUC- NO. OF ROUS IN THE DIM. STAIFMFNT FOR CjPHlrUHfKKiEX 
IN THE CALLING PROGRAM 

C= THE DAMPING MATRIX I 

PHI= THE MODAL MATRIX-EIGENVECTORS NORMALIZED TO UNIT MASS. 

UN - VECTOR OF UNDAMPED EIGENVALUFIS . 

ERf EI= VECTORS CONTAINING THE REAL AND I MAG . PARIS UF THE EIGEN* 
OF THE DAMPED SYSTEM. 

NTITLE= A VECT0R(DIM--80> CONTAINING THE TITLE OF THE RUN. j 

IPRT = A FLAG TO ALLOW PRINTING THE DYNAMICAL MAT . ( 1 - PRIN f , 0 NO 
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IV. OTHER SUBROUTINES USED 


ASMD2 calls NASA Library routine RESV to compute the eigenvalues 
of the dynamical matrix [A] . The computed complex eigenvalues are 
stored in array ER(I) and EI(I), which are the real and imaginary 
part of the Ith eigenvalue of matrix [A] . 
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V. LISTING THE SUBROUTINE ASMD2 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

t. 

f 

r 

c. 

(, 

L 

C 

I. 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


100 


200 

c 

c 


13 

12 


SUBROUTINE ASMD2 (N »NM» NRGVIC r C»PH I » UN » ER» El » NT1TLE* I PRT ) 


THIS PROGRAM USES THE ASSUMED MODE Mt 1 HOD TO COMPUTE 
COMPLEX EIGENVALUES OF A DAMPED SYS 1 EM. 

A FULL DAMPING MATRIX MAY BE USED ♦ 

THE PROGRAM IS A MODIFICATION OF PROGRAM CEIG 4/13/82. 


FORMAL PARAMETERS: 

N- N.D.O.F. IN I HE SYSTEM 

MM- NO. MODES USED IN THE ASSUMED MODE CALC. 

NROtIC- NO. OF ROWS IN IHE DIM. ST AT F ril'.NT FOR C f PH I , UK' , k R » F. I 
IN THE CALLING PROGRAM 
L- 1 HE DAMPING MATRIX 

f HI- 1 HE MODAL MATRIX-E1GENVEC IOKS NORMALIZED TO UNIT NASS. 

UN = VECTOR OF UNDAMPED EIGENVALUES. 

FR»Ll- VECTORS CONTAINING IHE REAL AND IMAG. PAR IS Oh THE F. L GENUAL. UE 
OF THE DAMPED SYSTEM. 

NT I f LF - A VECTOR ( DIM* 80) CONTAINING IHE TITLE OF THE RUN. 
i !•:*• 1 - A FLAG HJ ALLOW PRINTING THE DYNAMICAL MAT . < 1 "PRINT t 0 NO PRINT 

!> 

HI ME NS 1 ON CC MRUUCf 1 ) » PH I ( NROUC » 1 ) » UN(NROUL ) »ER(NROWC) * 

LI (NROUL) » N l 1 1 LE ( 80 ) 

D l Ml. NS ION DYEXAG ( 99> 99 ) , V(99»S) » UK (270) > A <99 >99) 


MAX 

NDOF -M 
N M l‘i • 2 >:> N M 


FORM DYNAMICAL MATRIX 


NN--2X NM 
DO 100 1=1 »NN 
DO 100 J= 1 » NN 
A ( I r J ) =0 . 0 

DYEXAC(I»J)=0.0 

DO 200 1=1 »NM 
IN=I+NM 

DYEXAC<I»IN)=-UN<I)*UN<I> 
DYEXAC ( IN » I ) = 1 . 0 
CONTINUE 


DO 12 J = 1 » NDOF 
DO 12 L= : l t NM 
DO 13 K=1 »NDOF 

A< J»L)=A< J,L)+C< J»K)*PHKK»L> 
CONTINUE 


56 


c 


DO 14 K=1 , NM 
DO 14 L=1,NM 
DO 15 J-l , NDOF 

15 DYEXAC ( K, L >=DYEXAC < K, L > -PHI <J,K)*A(J,L> 

14 CONTINUE 

C 

PRINT*, ********************* OUTPUT *##**###*#**##****• 
PRINT 11 

PRINT 10 » (NT ITLE ( I > , 1 = 1 , 80 ) 

10 FORMAT ( 80A1 ) 

? PRINT 11 

11 F URMAT ( // ) 

IF* IPRT.EQ.O) GO TO 210 

i. 

PRINT 201 

201 FORMAT( 1H1 ,///»* DYNAMICAL MATRIX IS */,/) 

DO 202 1=1, NN 

202 PRINT 205, ( DYEXAC ( 1 , J > , J = 1 , Ni! > 

205 FORMA! </»5X»5E15. 6) 

210 CALL RESVtHAX, NMD, DYEXAC, F.R«Ei, 0,0,0, UK, I ERR) 

C 

PRINT 1020, NM 

PR IN I*, 1 THE DAMPING MATRIX IS.** 

DO 666 1=1, NDOF 

666 - PRINT*, (C(I ,J>, J-l, NDOF) 

C 

;«)!(> FORMAT </» 10X, * ASSUMED MODE METHOD, USE *»I5»»* fUml •_.*/ '< 

PRINT LI 

PRINT*, 'THE DAMPED EIGENOALUk S t * 

<; 

DO 6J6 I 1=1, NMD 

JS3 PRIN 1 * , 1 1 , ER ( 1 1 ) , • * * , El ( XI ) * ’ J* 

PR IN I 11 
PRINT 11 

•.'O?? CUNUNUr 


l 
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VI. CONTROL CARDS AND TEST PROGRAM 

The test program TESTAS4 can be run through the use of the job 
file of Table 1. Table 2 is a listing of the test program TESTAS4 . 
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Table 1 Control Cards and Input Data for Sample Problem 


/JOB 1293B 

BPU » CM300000 » T1000 . 

USER' » 697121E . 

CHARGE»102428»LRC. 

DELIVER. 1293B D UVA 
ATTACH»FTNMLIB/UN=LIBRARY. 

GET f TESTAS4 . 

HAP ( OFF ) . 

FTN»I=TESTAS4»L=0tQPT=l»B=F. 

LDSET »LIB = FTNMHB. 

F. 

REWINDtOUTPUT. 

COP YE I » OUTPUT »0TTAS4 . 

REPLACE »0TTAS4. 

DAYFILE»L=DATAS4. 

REPLACE! DATAS4. 

EXIT. 

REWIND » OUTPUT . 

COPYEI t OUTPUT » X0TTAS4 . 

REPLACE » X011 AS4* - 
DAYFILE»L=XDATAS4. 


PACK » XDAT AS4 . 
REPLACE»XDArAS4. 

/EOR 

TEST: 2-DOF free system 4 ^ 


n , n 

*• » «■ 


1 . 3 »-. 5 » 0 » 0 < 4 - 

0 ^ 

0 » 1 ^ 


0.57 7» - 0.577^_ 
]. 2247>l4- 


h flz 

% DOF, I/J'cl 


■w, 

-f, 


0.81&»-0. 408 < 

0 

EOI ENCOUNTERED. 


U) k - &K‘ I 


-t 


/ 
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UVA 





Table 2 Listing of program TESTAS4 


C 

c 


c 

c 


p* 

PROGRAM TESTAS4C INPUT » OUTPUT ) 


DIMENSION G<88)»C( 99 ,99), F‘H 1(99,99), WN(99)»ER(99),EI<99> 
DIMENSION NTI TLEC80) 


1 

C 

C 

C 


10 

599 

C 


20 


C 

C 


5010 


99 


' 999 

C 

c 


c 

c 

c 

c 


DO 1 1=1 >50 

DO 1 J : =l » 50 
C(I»J)=0,0 

INPUT SECTION 

READ 10, (NTnLE(I)»I = l#80> 
PRINT 10, (NTITLE(I) ,1=1,80) 
FORMAT ( 80A1 ) 

READ *,N,NM 
PRINT *rN, NM 

DO 20 1=1, N 

READ* » ( C ( I » J > , J=1 , N ) 
PRINT*, (C<I, J) , J=1,N> 

READ *,IPRT 
PRINT * , IPRT 


DO 5010 I = 1 , NM 
READ*, UNCI) ,0<I) 

PRINT* , UNC I ) , G ( I ) 

READ*, ( PHI ( J , I ) » J = 1 , N ) 
PRINT*, <PHI(J»I),J=1,N) 
CONTINUE 

PRINT*, * N WN G* 

DO 99 1=1, N 
PRINT*, I, UNCI), G(I) 
PRINT*, -PHI' 

DO 999 J-l rN 

PRINT*, (PHIC J,I) ,I=1,N) 


CALL ASMD2<N,NM,99,C,PHI,WN,ER,El,NTnLE,IPRT) 
PRINT*, •1=REPEAT,0=ST0P* 

READ* , ITEST 

IF C ITEST » EQ « 1 ) GO TO 599 


END 
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Fig. A2.1 Actively Controlled 2-dof system 


61 




VI I SAMPLE PROBLEM 


Figure (A2.1) shows an actively controlled 2-dof free-free 
system. The contorl law is 


= -1.3 x^ + 0.5 

The input data to test program TESTAS4 is shown in Table 1. Table 
3 shows the output. 


Table 3 Output of Sample Problem 


******************** OUTPUT 


test: 2-DOF free system 


ASSUMED MODE METHOD* USE 

THE DAMPING MATRIX IS: 

1.3 -.5 

0 . 0 . 


THE DAMPED EIGENVALUES? 

1 0.+0.J 

2 -.3 402281689367+0. J 

3 -.4790459155316+. 971771301425 J 

4 4790459 1553 16+-. 97177 1301 425 J 


i: i t .j: :j: 


2 MODES 


Appendix A- 3 


Subroutine SPSTGN Documentation 


Damped Eigenvalues Using Reanalysis 


I. Introduction 

SPSTGN is a FORTRAN subroutine that can be used to compute the 
damped eigenvalues when damping is introduced into an originally 
undamped system through the introduction of a diagonal viscous 
damping matrix. Based on the modal data of the undamped structure, 
the damped eigenvalues are computed using an efficient reanalysis 
formulation. It is the purpose of this document to summarize the 
formulation and usage of the subroutine SPSTGN. 

II. Summary of the Reanalysis Formulation of the Damped 
Eigenvalue Problem 

Given an undamped system described by the following 
normalized modal data 


w i' i ~ 1 to L < Number DOF = NELEM 

where 

= natural frequency of the ith mode, rad/sec 
= OMEGA(I) 

(p. } = ith mode shape 

p . ^ = mode shape coefficient, jth component of the ith mode 
3 ’ = PHI(J,I) 


The normalization requirements are 

T r 0 if i ¥ i 1 

{p [M] {/>.} = < , “ 1 T ? f 

' i J 1 3 ' j il if i = }j 

where fM] is the system mass matrix. 

With the introduction of n damper at dof J ,J , . . .J 

c 1 2 n 


damped eigenvalue problem can be formulated as the solution 
characteristic equation: 


the 

of the 


detfl + £(s)6) = 0 


( A3 . 1 ) 


where 


(£] 


n x n condensed receptance matrix 
c c 
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with elements 



■Vi 


for i = 1 , n 

c 

for j - 1, n c 


The element R^_. are 


computing using 


R 


e fi 

L — 14 5-4— 

E 

£ = 1 2 2 
s + U) 

l 


(A3 .2) 


The following nomenclature is helpful in understanding the program 

L = no. of modes used in the receptance calculation 
( < number of DOF ' s ) 
n c = NSTAT = no. of dampers introduced 

= ISTATNO( i ) = ith damper location (DOF in the model) 
for i = 1 to NSTAT 


C 1 



0 


0 



condensed diagonal damping matrix 


c. = CAKRAY( 1 ) = ith damper value, for i = 1 to n 
1 [located at CT (or ISTATNO(i)] 

Two types of variation of damping constants are allowed . For no 

variation, set NCD = 1, IDOPT = 0, and CDINC = 0. 

A. IDOPT = 0, all dampers vary at the same time. For this case 
the initial damper values are c . = CARRAY( I ) . Then all the 
c. f s are increased by CDINC after each pass through the code. 
Tliere are a total of NCD increments as specified by the user . 

B. IDOPT = 1, only one damper value changes, all other values 
remain constant . The initial damper values are input through 
the array CARRAY( I ) , and the initial damper values are 

c . = CARRAY( i ) for i = 1 to n . 

1 c 

C^ will increase by CDINC each time the damper constants 

d are updated . 

C = CARRAY( n . ) + CDINC 
n , d 
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r» d = NODEPL = the number of the particular damper in 
ISTATNO(I) being varied. 

Ill Input Parameters 

111 . 1 Formal Parameters 

Description of the 11 formal parameters: 

NELEM: number of DOF’s in the original system 

NMOD : the number of modes requested for the printout < L 

RUNNO: your run number identifier 

IDOPT: o = all C values vary 

1 = only one value varies ( see NODEPL ) 

EPS : convergence parameter ( used in the test for convergence in 

the Muller’s Method root finding routine} 

NDAM : number of modified system modes desired ( < L ) 

NUMITER: maximum number of iterations per eigenvalue 
CDINIT : initial damper value ( units should be consistent with 

mass and stiffness, e.g. lb sec/in) 

CD INC : increment in damper vlaue 

NCD : total number of damper increments 

NODEPL: the number of the varied damper in ISTATNO( I ) . NODEPL is an 

element of 1,2, . . . ,NSTAT. 

Note: if IDOPT = 0. 

C = C iQ + Ac for L = 1,2 NSTAT 

if IDOPT * 1 

c. » c. i # NODEPL 

C. = c 1 + Ac, i = NODEPL 
l 10 

111 . 2 Common Blocks 

. In addition to these variables, the following variables enter 
the subroutine through the common blocks and must be evaluated 
prior to calling SPSTGN. Also, the common block set contained in 
the file COMSPST must be placed in the main program. 

NTITLE( I } : a holerith array with a title of the job 

L: the number of modes used in the receptance calculation 

OMEGAf I ) : a vector containing the L values of the original w ^ 

PHI(J,I): a matrix whose L columns are the modeshapes of the 

original system (DIM: NELEM x L) 

Note: These modeshapes must be normalized to unit 

generalized mass. 

NSTAT: number of stations where dampers are placed 

ISTATNO( I ) : an array giving the location ( DOF ) of the attached 
dampers I = 1 to NSTAT 

CARRAY( I ) : an array containing the damper value for position 

ISTATNO(I), I = 1 to NSTAT 
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IV. Using Subroutine SPTGN 


The data required are input through formal parameters ( Section 
III.l) and through the common block set COMSPST, into SPSTGN (Section 
III .21. It should be noted that mode shapes have to be normalized 
to unit generalized mass. That is, 

{ 1 P i > T [M]{P i } = 1.0 


V. SAMPLE PROBLEM 


The Model is the NASA 88 D.C.F. grid. 

The model data is read from TAPE 12 is NASSAGG2, which contains data 
for all 88 modes, including 3 rigid body modes. 

Four dampers are placed at D.O.F. 1, 11, 77, 88. 

The initial damper values are 0.2. 

The approximate receptance elements are to be calculated on the 
basis of using L( = 30**) modes of the 88 modes. 

The first 2 modes are to be printed out ( NMOD = 2 ) . 


All the damper values are to vary ( IDOPT = 0 ) with increments of 
CDINC =0.1 for NCD 5 trials. 

The convergence criterion for eigenvalues is ldet(A)l < € = 10 


The test program TESTSP is designed to interactively accept the 
above data and then call subroutine SPSTGN to do the 
calculations and print out the results. 


The following input values were given ( the remark number 

corresponds to the number in the following TESTSP printout ) . 


Remark 3 TEST SPST (title ) 

2 88, 2, 30, 1 (88 DOF, 2 modes to be printed, use 

30 modes , Run No . = 1 ) 

3 0 ( all damper values change ) 

4 The modal data was read from TAPE 12 NA.|AGG2 

5 4, l.E-S, 2, 50 (4 dampers € (EPS) = 10 . 

computes 2 damped eigenvalues, max 50 iterations) 

6 1, 11, 77, 88 (damper locations) 

7 0.2, 0.2, 0.2, 0.2 (initial damper values) 

8 0.2, 0.1, 5, 1 (initial damper values, increment, 

5 increments, all dampers vary) 


**An arithmatic underflow condition will appear if a larger L is 
used, we have had such an experience by using L = 88. This is 
due to the formulation used in teh receptance calculation. It 
can be improved if so desired. Our experience indicates that 
using L = 30 provides very good results. 
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VI. Program Listing and Sample Run 


The subroutine SPSTGN is listed in Table 1. Table 2 lists a 
sample main program TESTSP which calls SPSTGN. The output of the 
sample problem is given in Table 3 . The control cards necessary to 
run this sample problem is shown in Table 4. 
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Table 1 Listing of Subroutine SPSTGN 


SUBROUTINE SPSTGN ( NELEM » NMOD .RUMNO » IPOPT » EPS t NPAM » NUMITER > 
CI’INIT f CD INC r NCD » NODEPL > 


p 

p 

r 


3PSTG WRITTEN BY DR * A. PALAZZOLO AND 
COMPUTER CODE TO GENERATE THE COMPLEX 
R E A N A L Y S 1 3 F 0 R M U L A T • 0 r! , 


DR. B , P, WANS MAY 7 1981 
EIGENVALUES USING 


C 

r 


C 

r 

p 

p 

p 


INPUT CARD DESCRIPTION 


nelem: no, 
nmodi no. 


OF ELEMENTS 

OF MODES REQUESTED FOR PRINTOUT 


AND DAMPER SENSITIVITY 


10 AND 0 L 


ANALYSIS AND PLOTTING 0 
EUNNO: RUN NUMBER IDENTIFIER 
IDOPTt 0= ALL C VALUES VARY 

1 " ONLY ONE VALUE VARIES 
E P S - C 0 N V E R G E MCE P A K A METE R 
HDAM: NUMBER OF MODIFIED SYSTEM MODES 
f ! U MITE & T M A X J M U fi M U M ft E R OF IT F R A T I C N 3 
AL I AM REF VALUE FUR ROOT 
E VILL OVER IDE T HP. VALUE 


CD I NITT [NT 


DF SIRED S L-2 
PER EIGENVALUE 
LOCUS PLOT LB 


AND Q2\ 


ccr 


THIS 


Hi 


/ T M 

IN POSITION NODEPL ? ENTERED 


INTO 


CARR AY 

CDINC: DAMPER INCREMENT FOR ROOT LOCUS PLOT LP.SEC./IN 
NCDJ TOTAL NUMBER OF DAMPER INCREMENTS FOR ROOT LOCUS PLOT 
NODEPL*, THE NUMBER OF THE VARIED DAMPER IN ISTATNO(I) 

I.E. NODEPL IS AN ELEMENT OF 1 ? 2 j , , » NST AT , 

THIS TELLS WHETHER NODEPL CS FIRST • SECOND ETC, IN IBTATNO 


N ADM": LON THE FOLLOW 
p r i m u n >i <( i o r : 


ENTER THfc' SUBROUTINE THROUGH THE 


■ 0 P P F. C E P T ANSI 0 A L C U L A T 1 0 i ' 


l ; ; in f pp D D.E : 

I i , I, i ! 0 ( J ■ : THF ! 1 S T A f DAMPER LGOATJUf!'- (NODE NUMBERS) 

L A R R A • J ; ’, THE DAMP!- R VALUES CORRESPONDING TO IS1ATNU' I *> 

NST A T: NUMBER OF STATIONS WHERE DAMPERS ARE ATTACHED 
RHJM.J,: AN ARRAY CONTAINING THE UNDAMPED SYSTEMS EIGENVECTORS 
NORMALIZED 10 UNIT GENERALIZED MASS 


0 M r - 
NT 7 


I T 


A VICTOR CONTAINING THE FRFQUF.NC. 


, i 


HA ARP AY CONTAINING THE TIT 


T Tl IT 


3 OF THE UNDAMPED SYSTEM 
OF THE PROJECT 


INTEGER lETATi-N- 

i ■ r m r: m •: r I > m i P on; i , \ r. , „ r i two.' i a a 

CfJMPLE< Cl GEN » SLAM IN » VLAM > RLAM t CUEIG » SLAMNXT f CMCOEIG » F I I 
rOMPLEX FH AT 
! f; 0 1 C A l TEST 
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Table 1 Continued 


C 

r 


COMMON /BLK 1 / AE J ( 50 ) f ASL J C 50 > ? AT J ( 50 ) » ARH J (50) » ARQ < 50 ) 
COMMON / BLR 2 / SK(1f1' 

CDMMCN/BLN3/SSTR(88) fSSROT(BS) 

C0MM0N/BLK6/UR (S3) » Z< 1 , 1 ) .PHI ( SB ?33 ) r OMEGA (83) ?DYN( 1 1 1 ) 
6 EM(lvl) 

C0MMQN/BLK7/VLAM (3) ? RLAM (3 > » COEIG (20 ) 

C0MMQN/BLK8/L 
CGMM0N/BLK10/SMINV( 1 r 1 ) 

COMMON/BLK 12/UREX ( 105 ) fUIEX( 105) 

COMMON/BLK 13/ CMC OEIG ( 105) 

COMMON/BLK 14/SENSIK 1 t 1 ) 


COMMON/ BLK1! 


i( 1 » 1 ) 


COMMON /BLKI&/PMR1 

COMMON / B L K 2 0 / IAUTOfREMAXf A I M A X 

COMMA N / B L K 2 1 / 3 T 0 M ( 3 3 ) 

COMMON/BLK 22/FACT 

COMMON / B L K 2 3 - ISTATHU (16) ? NST AT » CARR AY ( 16) 
nQMMQN/HLi;2'l 'R1VEC(16> 

COMMON /BLK25/IME 
C 0 MM Of! / BLi; 2 6 /FHATUOrl 0 ) 

COMMON / B L K 2 7 / NTITLE ( 30 ) 
t f K t Y Y Y t 

U! Y Y Y % t Y IT 1 t'YY n. * 

CM =1 v 0 

P L T A - 1 ♦ 0 
I D I A 0 - 0 

READ t , NELEM j NMOD f NRIG 7 L f I FLAG » R f BETA 7 RUNNO 
PRINT 1 f NELEM f NMOD f NRI G f L f IFLAG f R f BETA f RUNNO 
FORMAT ( 51 10 f3F1 0.3) 

N= NELEM 
lW-NELEM+1 
READ Y f J f!OPT * 

PRINT ?5f[N 0!T 
TE' JllOPr.EQ.C-; GO 10 MO 
GO ■ 0 250 
i m f i r j m ij t 


BOP I 

■ non 
go ro 


’ . i ; •_ ; > \ 

!■’ !- A H ( ! 2 7 0 1 I ; 0 M E G A < I • r G M 
5 COM ( I ) OMEGA C I ; 

C READ ( 12 7 '“01 i ) •' F’H I ( J 7 I ) » .} • 1 7 i! ) 

501'' CONTINUE 
5011 FORMAT ( CE16 .-') 

5' 11 FORMAT ( 5 V ; 2 E 1 6 * 9 '■ 

50 ! 0 210 

'50 l Of! I T fi’JL 

DO 5 [ 1 ? : I • 1 f ! 

C P E A I' + f S S T R ( I > f 5 3 R 0 T • I 7 

'5 T l- HIT s T R ( I ) r £ 3 R G T •' I ) 

C 6 FORMAT C2F 10. 3) 

C 5 CONTINUE 

r HO 10 I 1 r NELEM 
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Table 1 Continued 


: READ * ? AE J (ID ARL J ( I D AT J ( I D AEH J (ID; 

PRINT 15 ? AE J (I)r ARL J (DrATJ(I)r ARH J ( I ! 
15 FORMAT' 5F 10*2) 

10 CONTINUE 
40 CONTINUE 

READ * r MST AT ? EPS i NDAM * NUMITER ? ID I AG 
P R I NT 2 0 y N S T A T - E P S 7 N D AM 7 NUMITER ? ID I AG 
20 F 0 R M A T ' 1 1 0 7 E 1 0 . 3 7 3 1 1 0 ) 

READ 1 1 5 7 T M D P L 1 


C 

C 


r 

u 

p 

r 


/ / 


j j 
L 


READ * . ( 1ST A TMO ( I ) 7 I- 1 ? NST AT ) 
PRINT 757 (ISTATNO(I) ;I-DNSTmT> 
FORMA T< 1615) 

READ * ? ( CARRAY ( D ? I 1 ? NSTAT > 
PRINT 7 7 7 ( CARRAY ( I '< rI--DNSTAT> 
FORMAT (3F 10.3) 

READ K 7 CDIMIT 7 CDINC , NOD 7 NODEPL 
P r j vi t ^ 5 7 C D I N I T C D I N C ’ H C D ? NO B E P 1 
FORMAT ( 2 F 1 0 , 7 ? 2 1 1 0 ) 

FORMA f ( 2 1 1 0 7 2 F 1 0 ♦ 3 > 


r 

r 

i - 

\u 

C 

r 


FORMAT'" 1 H - ? 1 0 X 7 DIME 1 - DK10.2) 

tfN ft ftt* 

K** "D ****** Y 
RE ANALYSIS SECTIOi! 
*************** 

C A R P A Y < NODEPL > = CDINI T-CD I N C 
CALL INPR2 ( L 7 EPS 7 NDAM ) 

DO. 3 0 D-D NOD 

c; A P K A Y t i ! 0 D E P L > • C A P R A Y (NO D E * L • -i C DIM C 
CAL! i UPP2-: L 7 \i?- - NDAM > 

DU 50 DDE TO i ; NDAM 
FACT -1 . 0 

CALL I M T T I AL •" S L AMI N 7 [ NL'F. ; 0 7 PET A • I ; 

I U i J TER" ! -NUMITER 
UAL! tl ’ 1 L. E P •' U L A M ? I '/ T > 

] m :: Ti ^ r p c j j* Mf ’X I I U D L I 0 • 

CA! L C Ci i 1 U C H I > ' [ TLR 7 LPC - ! LSI ) 

I M I-- T r p ■ ■ 

1 F ( 1 L T ; 00 10 70 
40 CONTINUE 

70 CALL EI5ST0R' INDEIG) 

CAL!. M EXMOE’F ( CL AMIN r IL'DEIG > 

COflTTNU 1 - 

'Yi! L APEISOT 'MPAM) 

CALL AFEIPFD NDAM; 


PR TNI 103 ? T 

IF(IPOPT.EQ.O) GO TO 1100 
CARRAY ( NODEPL ) =C ARRAY ( NODEPL ) +CDINC 
CO TO 20 
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Table 1 Continued 


10 *' DO 1000 I B = 1 r NSTAT 

0 0 0 C ARRAY (ID >=C ARRAY ( ID)-! CD INC 

103 FORMAT ( 1 H » 2 X»*TIME 3 " t,F 10 . 3 ) 

30 CONTINUE 

RETURN 

END 


SUBROUTINE INPR2 (L; EPS » N D A M ) 

INTEGER ISTATNO 

COMMON/ BLK23/ ISTATNO (1 6 ) 7 NSTAT ? C ARRAY (16) 

PRINT 1 

FORMAT < 1 H 740X7* BEGIN REANALYSIS PORTION OF CODE* > 
PRINT 20 

FORMAT ( 1 H 0 7 5X ? * N 0 D E * r 1 OX 7 * DAMPER VALUE * > 

DO 30 I - 1 f NST AT 


PRINT 


TAT N 0 ( I ) : CARPAY ( I > 


3 1 '. FORMAT ( 1 HO? 6 X ? 12 ? 1 5 X r E 1 0 » 2 ' 

50 CONTINUE' 

PRINT 5 y L r E P 8 

3 FQRfiAT'lH ? 10 X"*MO, OF MODES FOR RECEP. CALC . ; L 

'7 y L 1 * w / 


* 




t;< m r r> r c i c c t c •z 'i 


PRINT 10’NDAM 

■ r ■ I c- !,| 

l I 1 • C I . 1 

END 

SUBROU TINE I N I T I AL ( SLAM I N - I NDE 15 t PET A ? I ' 

CALCULATE THREE STARTING VALUES FOR MULLER ITERATION 
INTEGER ISTATNO 

COMPLEX C I GEN y SL AM I N 7 VL AM ? PLAN 7 COE I G j SL AM ? RL ■ CMPLX t EE 
COMMON/ BL i* 7/ VL AM ( 3 ) 7 RLAM ■ 3 > 7 C 0 £ 1 0 • 20 ) 

COMMQN/DL! S/!. 

•COMMON 'EM MCI 'ETON (38 • 

1'ijtiw ft>,i 'Tii i, 22/ FACT 

; , 1 ‘ r1 M Or 1 ' L ! . ’ 3 / IS" AT NO ■' i 3 > 

COMMON 'DL L20 ' I ML 
MPTG n 
i R ] AO 0 


7 C A R R A ' 


i r c i ; ! r I ir r 


; M, L ' 


if, ! 


r r ! ni r * 


r.ti ip i 


! F ■' E " i AO , 2 C . 0 
I M< L NT 17ICL 

FORMAT* 1H 7 4 0 X 7 {'ITERATION FOR EIGENVALUE S ? I 0 > 
CON TI M UI". 

CZ-STi’MClCC: 

1 " r ;< p j y ' - 0 . C r 7 E ’■ 

1 1 • r r .* v r 7 7, 
r • ■ r . •- n 1 1 am r m - c •- 

I • J. <• I. V 4. J ■... I I 1 I l O. I \ I I _ 


j J 

i r 

i!E 

. ! GLAM 

I M ~ C 0 E 

IG( 

I N D E I G 

V ] 

A M < 

- '1 

--SLAM IN * 

( 1 . 0 7 " 

, 00 

1 ) 

V L 

am •: 

2 ) 

- 1 » 001 *S 

LAM IN 



VL 

AM ( 

1 ) 

-SLAM I M 




0 1 

A M = 

VL 

AM ( 3 ) 




C :*i 

LL 

POFLACSLAM 

? RL 7 IN 

DEI 

G) 
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Table 1 Continued 


RLAM ( 3 > ~RL 
SLAM-VLAM ( 2 ) 

CALL ROFLA ( SLAM ? RL * INDEIQ ) 

RLAM ( 2 ) - RL 
S L A M - V L A M ( 1 ) 

CALL ROFLA ( SLAM » RL » INDEIG) 

RLAM ’ 1 )-RL 
AB-C ABS ( RLAM i 1 ) ) 

FP-L 

TP- 1 ,0/TP 

cb=ab**tp 

F ACT- 1 , 0/CB 

RLAM t 3 ) =RLAM ( 3 > / AB 

RLAM t 2) -RLAM ( 2 ) /AB 

RL AM ( 1 ) - RLAM ( 1 > /AB 

RETURN 

END 

SUBROUTINE M U L L E R ( S L A M N X T ) 

CALCULATE NEXT EIGENVALUE GUEFS 


B V M U L L E R 

i n 


COMFL.EX. SLAMNXTy VLAM y ELAM? CQEI Gy CSQRT 7 FZJ3 
6 F J 1 J 2 J 3 y U y 3 Q y A C 1 r A C 2 y D E N 
C0MM0N/BLK7/VLAM ( 3 > » RLAM < 3 ) ? CGEIG ( 20 ) 

F 7. J 3 .12 = ' RLAM ( 2 ) -RL AM ( 3 > ) / < VLAM ( 2 > -VLAM ( 2 > • 
F 2 J1 J2~ • RLAM < 2 ' -RLAM ( 1 > ) / ( VLAM < 2 ) -VLAM • 1 < ) 
F 2 J 1 J2" < RLAM ( 3 > -FLAM (!))/( VL AM ( 3 ) -VLAM ' 1 > 

F J 1 J2 J3 • ' FZ J3J2-F Z J1 J2 ) / < VLAM ( 3 > -VLAM ( 1 > > 
W=FZJ1 J2+FZJ1 J3-FZJ3J2 
SQ = im2--4 . 0*RLAM< 1 }*FJ1J2J3 
SQ-CSnPT(SQ) 

AC 1 -WtS i! 

A C 2 = U - S 0 
rt- CABS < AC l ) 
r':-CABS t AC'"/ 

•: r • n , (- e , r ? • 

I F ■' r 2 . L' I , I 1 ; 

SLAMNH VLAM 
T ETURrl 
END 

' ; 1 J B R 0 ' 1 flilF '."-'‘DATE ( 3LAM.N/.T 


METH03 
FIJI j: 


^ 7 


DEN- AC 1 
DEN -AC 2 
1 ) -2 -01 FLAM ( 


• f 1 m i l- 


M ' 




LClMMON/FL R7/VLAM • 3 ) - RL 


: n 1 


"PDA T E THE LAMBDA ( VLA 
M 3 ; - VLAM ( 2 ; 

M- 2) = VLAM< 1 ) 


/ n r ' l 


- j- ; d T T p, \ 

•: (J E J G - h ! 
> y CUR [ G ■ 
> AND i-‘ 


20 

OF 


iMBDi 


ARP AY 2 


: ; i / 

VLi 

"I AM ( 1 > -GLAMNXT 
C A!. L K 0 F !. A ' 3 LA 31 NX I y RL t I N D E I G ) 

RLAM -'3/ PI A M ’ 2 '• 

p! t ‘j r ; . 

RLAM ' ; • - RL 
RETURN 

l: (id 

SUBROUTINE CONVCHK ( I TER > EPS , TEST ) 

LOGICAL TEST 

COMPLEX V L A M y R L A M j C Q E I 0 j S L A M * TEST 1 r I FST2 ; TE' 
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Table 1 Continued 


C0MM0N/BLK7 ' V L A M ( 3 ) ? RL AM ( 3 ) f CQEZG(20> 

COMMON ''BLK25/IME 

TEST FOR CONVERGENCE IN MULLER METHOD 
I D I A G = 0 
EPS2 -0 ♦ OILERS 
SL AM-VLAM ( 1) 

T E C T 1 ~ 0 L A M ( 1 ) - V L A M ( 2 ) 

TES I 2~RLAM ( 3 • 

7 L S 1 3 ~ V L A M ( 3 - V L A M ( i > 

CHECKl -CABS (TEST! > 

CHECKS -'CABS (TEST2) 

CHECK 3~CABS(TEST3) 

TEST - .FALSE* 

IF ( CHECKl . LT . EPS . AND . CHECKS * LT * EPS ) TEST- , TRUE . 

IF ( CHECKl . LT . EPS 2 ) TEST= . TRUE * 

IF < CHECKS, EQ, 0,0) TEST~=, TRUE, 

IF ( I D I A G , E 0 ♦ 0 ) GO 70 5 

PF; I NT 1? LTER? VLAM< 1 ) 7 RLArl ( 1 ) ? CHECK 1 - CHECKS ? CHECKS 
l FORM A H 1 H ? 1 X r 1 2 r L X ? E 1 2 < 5 p 1 X r E 1 2 , 5 p 3X ' E 1 2 ♦ 5 y 1 X ? E 1 2 . 5 ? 3 ( 5X r E 1 2 . 5 > ) 

5 CONTINUE 
RETURN 
END 

3UBR0L 1 * * * S TINE EIG3T0R( INDEIG ) 

C 0 M P L E X V L A ii ? P ! . A M ? C 0 E 1 0 

S 1 0 7 ■: E M 0 1 1 1 F [ E D 0 Y S 7' E M EIGEN v A L U E 
C 0 i 7 M U i ’ / B L K 7 / LAM ( 3 > ? RL An ( 3 > r COEIG ( 20 > 

COEIG( INDEIG) =VLAM< 1 > 

RETURN 

END 

SUBROUTINE NE X MODE ( SLAM I ?! 7 I NDE I G ) 

C 0 M P L E X S 1. A rl I N p V L A M ? F L A M » COEIG 

COMMON / B L IK / V L A M < 3 ) 7 RLAM ( 3 > ? COE I G 

C Q M M 0 N / B L K 2 C / I ! ! 7: 

THIS 3UPP0UTINF. IS C U F’ RENT!, f NOT BEING USED, 

SUBROUTINE INTTFAL IS NOW MAKIiiG THE INITIAL GUESS 
ICE INDEIG 

i D 7 A f ’ -0 
N!' ! 0 =0 

r f / ; 1 1.- [ - , i' r i . I j I r F - T 7 1 D E I G -!• J 

LALLULATL INITIAL CENTRAL EIGENVALUE GUESS TCR MX 7 MODE 
'SLAM r H" 1 ,00 1 TCOi: 3 G ( I NDE I G > 

PRINT 0 7 ICE? L ME 

C FORFiATdH 7 0 0 X ? T E I 0 E : 1 V A LUE K » I 3 7 G X * TNG. OF ITERATIONS - ¥ - 1 5 ■ 

RL f 1 ’• R J 1 

S i ID 

SUCROIJ! [HE ROFI.A 'SLAM-P! -JMUEIO; 
rc.-.cc r = T AT flf 1 


/ r i n 


CCXPLE . 


) !_ , : r 


PLAn ’ LOEIG 7 FI I ’ 7 ION. 


. s n ' 


SUIFOUI IrlE T 0 EVALUATE IMF ITERATED FUNCTION 
COMMON / tLl S 'UR •' 3 3 ) - 7. ( 1 ;1), PH 1(38-88)7 OMEGA ( 33 ) ? DYN < 1 - 1 ) i 
C 0 M M 0 N / B L K 6 / U R ( 8 3 ) 7 Z • 1 7 1 ) 7 PH I ( 88 7 S3 ) 7 OMEGA ( SB > » D Y N ( 1 7 1 ) 7 
6 SM( 1 7 1 ) 

COMMON / B L K 7 / V L A M < 3 > 7 R L A M ( 3 ) 7 COEIG (20) 

C0MM0N/FLK3/L 
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Table 1 Continued 




COMMON /RLK22/F ACT 

COMMON /3LK23/ 1 STATNO ( 16) ? NST AT » CARRAY (16) 
COMMON /BLK26/FH AT ( 10? 10) 



C 

PRINT 220? FACT 



V 

FORMAT*; 10Xr*FACT = *f2El6.5) 



c 

PRINT 100? SLAM 



100 

FORMATdOX?* SLAM + *?2E16.5/> 




NRIG-0 
BETA- 1 * 0 




CALL RECEP(SLAM) 



c 

PRINT 100 "SLAM 




HO 20 1 = 1 ? N S T A T 
DO 30 J = 1 ? N S T A T 



20 

BBC ( I ? J ) = F H A T ( I ? J ) % C A R R A Y < I ) *SLAH*BETA 
BBC ( I ? I ) =BDC ( I ? I ) + C 1 . 0 ? 0 * 0 ) 



00 

CONTINUE 
1ST A f -NST AT 

IF < NST AT . Eli ♦ 1 ) G = BBC < 1 ? 1 • 
IF ( NS F AT ♦ ER , l ) 00 TO 00 
C; A L L CDETERtlST A T ? BDC ? 0 ) 



’JO 

CUNT IN UP 



p 

PRINT 100? SLAM 



r 

PPINT 310*0 



5 ! 0 

FORMAT ( 1 0 X * Y 0 !• ? 2E 1 6 , 5/ > 

p i r 




I T l_ •-* 

IF (pJRIG.EQ. 1 ) RL-G/3LAM 



C 

PRINT 3 1 0 ? G 



c 

PRINT 330 ? BETA 



330 

FORMAT ( 10X?*BETA=*»2E16.5> 




DO 1 I = 1 7 L 



c 

PRINT 3 4 0 ? I 


/ 

7 VO 

FOF MAT ( 1 OX ? T: Id? IS) 


1 

j 

RL -P!. ) ( SLAMttCT < BE T A ft 2 ) +OMEG A • I '• «2 : tF ACT 
IF' E M D t 1 G , E 0 , L '■ 00 TO 10 
[NT • TiJDF.IG-1 
00 2 [ -1 - LNT 



f ' 

PFTNT 200 ? I ? f'OE I G < i ' 



I. {) 

P OkiiAT • 10/ ? U ? COL 10 (DMT; 2L16 . C ) 




;:l '• 1 1 ' f CLAM-onr lm [ ;• > 
: ■’ A ,, DiAG* COE EG*' E • ■ 


1 


El- ■' KA-EO.O.O - GO TO 2 



( 

PR] NT 1 1 0 ? COM J 0 ' C ,_I E I G ( I ,* ) 



1 1 0 

F 0 R M A T < 1 0 X 7 6 COM J G , 2E 1 d . 3 ) 




RL -PL,' •; SLAM-CON JG < CGEIG ( I ) > ) 




CONTINUE 



l A 

continue 

r-:L f ! ! F ; M 




s y.u 

0 U k F U U ! E M E E C E P * S L A M ) 



L 

G 0 M F U T F P E C F. P TANG E E ; 



L 

MODFC NOF;MAL I ZED TO UNIT MODAL MASS 




INTEGER ISTATNO 
COMPLEX SLAM ? F 1 1 




COMPLEX D E N ? F H A T 




i 

i 
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Table 1 Continued 


C 0 M M 0 N X B L K 6 / U R ( S8 ) ? Z ( 1 ? 1 ; . PH I < 88 >88 > ■ OMEGA • 38 ' . DYN (1:1 
C 0 M M 0 M / B L!\ 6 / U R ( 8 8 ) »2( 1 >1 ) ’PHI (88: 83) ' OMEGA ( 88 '* ? DYN •’ 1 > 1 
6 SM ( 1 ? 1 ) 

COMMON XBLK8/L 
COMMON /RLK16/PMR1 

1 OMhON/PL!' 03 7 1ST A TNG 16; -NET A T ? C AREA Y ( 16 ' 

COMMON /BLK26/FHAT ( 10? 10) 


DO 5 I - 1 ? N 8 T A T 
TP- I STATNOi l ) 

DO 10 J ~ 1 ? N 3 T A T 
P HA T ( I ’ J > - ( 0 * 0 > 0 ♦ 0 ) 

JP - I. C 'T ATNO ( J N 
DO 20 K - 1 > L 

20 FH AT ( I : J ) -FHAT ( I > J ) f PH 1 < IP u ) fFHI < JP ’ L ) * 

6 'SLAM**2*(BETA**2>+0MEG A ( K M' * 2 ) 

10 CONTINUE 
C CONTINUE 
NET URN 
F ND 

SUBROUTINE APE 1 EOT (rlDAN) 

C C M P L E X V LAM? P L A M > C Q E I G ? D E : ! 

C (5M H 0 N ' D L I > 7 / U L A ri < Z ) ? R L A • 3 ) - C 0 E I G ( 2 0 ' 

A P RANH F ,1 0 D ! F I ED 3 Y 8 T E M E 1 0 E N V A LUES (A; P P F 0 , ( ' 
A C 0 0 KLI.LN3 10 INC R E A 8 1 : ! G I M A G , P A P T = , D E - T I M E 
PET A = 1 t 0 
DO 10 I- 1 >NDAM 
10 COEIG(I) -COEIQ(I)*BETA 
NPC-NDAM-i 
DO 1 I-lfNPC 
. 1 8 T 0 R - 0 

t ft f ?{ - 1 00 0000 , A 

! 3 T - 1 

m : j- TOT /NAM 
AT A EUAGCCQEibC J) '■ 


; i ! h y \ , \ 

i : . *, i i l 

I I r I I - • 1 


t.TMFM) TMTNAf 


! i- (AT . ! !i, HUM) 23T0i-' = J 

i*' ."tfi i '.nil. 

. - t < i . i 

DEN -I'Ni- £0(1' 

CNF I G • i ) LOEIC ( .iSTGR ■ 
i, HE LG ( J3 TOR ; -DEN 
1 ‘'ONI Ir n !L 

i,-U' T It p / 1 

!- ND 

VJt POD : ! NE APEIPRT (NDAM) 

CUhi'L’E •; AM ? PLAN, - COE I G 

L f, MMON/M ! 7/UL AM ( 3 ) - REAM ( 2 ’ COE 1 3 ( 2C 

P R [NT T H c M 0 B I P 1 E D 3 Y TEN: A P P P 0 X . EIGEN V A L U E 0 
-PINT L 

[ FORMAT ''IN ? 22X 7 * MODIFIED SYSTEM APPROXIMATE EIGENVALUE: 
6 *ZETA* ? 13X?*0MEGAN SEC-lt) 

H K I G ~ 0 

DO 10 I - 1 » M D A M 
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Table 1 Continued 


0 

C 


IF<r!RlG.EQ.i> ICE=I+i 
OMEOAN -CABS ( CQEIG ( I ) > 

Z.£ T r*i--RtftL ( COEI 3 •' I ) ) /0 ME GAN 
PRINT 20 « ICE > COEIG ( I > r ZETA i QMEGAN 
20 FORMAT < 1 H » 1 OX , 1 5 > 1 OX f E12 * 5 ? 3X j E12 , 5 ? 1 3X » F 13 . 4 r 1 IX » E 12 . 3 ) 
10 CONTINUE 
RETURN 
END 

SUBROUTINE CDETER ( NST » A ? D > 

DETERMINANT OF AN NST BY NST COMPLEX MATRIX WHERE NSTM 
COMPLEX A(lOflO) i D 
COMPLEX PMULT»BIQ»DET»BS»TEMP 
NS 1 2 = NST 
I COUNT 0 
NIJMSYS-NSTZ-1 
DO 14 I - 1 « NUMSYS 

SCAN FIRST COLUMN OF CURRENT SYSTEM FOR LARGEST EL EMMT 
CALL THE POW CONTAINING THIS ELEMENT ? ROW NBGRW 

HN 1 ~ 1 

BIG" A •' r • • > 

Mf(RDi,i r 

< A* >/ I « «. 

DO U J NNfNSIZ 

IF (CABS ( BIG ) -CABS ( A< J» I ) ) ) a ? 5 7 G 
6 BIO-A(Jfl) 

NBGRW- J 
5 CONTINUE 
BG ~ t .O/BIG 

• SWAP ROW I WITH POW NT GPU UNLESS I -NBGRW 


IF- NBGRW. NF- , ■ ; [ CQUN T “ICOUNTT 1 

IF -'NBGRW- IT 7»10» / 

I WAR A MATRIX ROWS 

DO -J J-UiGI' 

T F M P - < 1 N B G P U i J - 
A ( N BG C -' W i J A ( L 7 J • 

A ( J * J ) -TEMP 


ELIMINATE UNKNOWNS FROM FIRST COLUMN OF CURRENT SYSTEM 
10 DO 13 K NN 7 NS I 7 

COMPUTE PIVOTAL MULTIPLIER 
PNULT =-A< K 7 1 ; *PG 

APPLY F'MULT TO ALL COLUMNS OF THE CURRENT A MARRIX ROW 
HO 11 J- NN 7 NS I Z 
L 1 A • !• ? J } -RHULT f A ( I , J ) -i A ( 1. 7 J ) 

13 CONTINUE 
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Table 1 Continued 


1 4 corn [f'JE 

E> = ( 1 . 0 <• 0 ♦ 0 j 
DO 100 1 = 1 j NS 1 2 
100 D=D*AfI*I> 

JC0UNT-( IC0UHT/2)K2 

IF ( JCOUNT « NE » ICOUNT ) E' = -D 

RETURN 

END 

SUBROUTINE HEADMG 
INTEGER ISTATNO 

C 0MM0N/BLK23/ I STATN9 ( 1 6 > r NST AT ? CARRA Y (IS) 

PRINT 1 

l FORMAT! 1H1 » 40X» + E I SPACK { EXACT ' SOLUTION OF THE MOD 
PRINT 20 

20 FORMAT ( J H 0 ? 5 X *> INODE? ? 1 0 X ? * 0 A M F' E R VALUE 1* 

DO 30 1 = 1? NST AT 
PRINT 35 ? ISTATNO*: I ) ? CARR AY < I - 
55 FORMA H 1 HO ? 6X ? 12 ■ 1 6X ? E10 . 3 ) 

30 CONTINUE 
RETURN 
END 

SUBROUTINE LABRT ( I 3U ? LHQL - INX) 

DIMENSION LH0LI5) 

LOGICAL PS*TS 

I F •: ( I SU ♦ ECU 0 > . OR . ( ISW ♦ GT . 5 ) '* RETURN 
GO TO ( 1 y 2 y 3 y 4 r 5 ) ISW 
DATA NP/ 10/ y PS/ . TRUE ./US/. FALSE . / 

IF (PS. AND. cNP.GT.O) ) PRINT 27-LHQLrlN/ 

FORMAT ( 1H0 t?'.‘ i 5A10 » 3X ?0i ) 

MP-fip~l 

1 F ( f *3 ; CALL EYl! 


;<r- f n p > i 

i ; p I it" , 

,\i ; . -1,1/ 

-E. : URN 

c ■; [".'IJL 

| i v 1 1 .* > : 

I Nl! 


F I ED 
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Table 2 Listing of Program TESTSP 


PROGRAM TESTSP ( INPUT ■> OUTPUT ? T ARE12 > 

SPSTG WRITTEN BY DP, A. PALAIZOLO AND DR * B , P, WAND MAY ? 1981 
COMPUTER CODE TO GENERA I E ROOT LOCUS PLOTS FOR A EULER BERNOULLI BEAM 
MODEL WITH MULTIPLE DAMPERS - 

THE BEAM IS SUPPORTED BY TRANSLATIONAL AND ROTATIONAL SPRINGS, 

THIS CODE WAS DEVELOPED UNDER THE NASA LANGLEY SPACE STRUCTURE CONTRA! 
INTEGER I2TATN0 
DIMENSION IROH(IOO) ? TEMP < 100 > 

C 0 M P !. E X C I G E H . S L A M I N » V L A M - R L AM.-CQEIGfS !.. A M N X T » C M C 0 E I C- ? F 1 1 


1 0 0 1 C A L T F T 

C lJ M M 0 N / £ L i\ 1 / A E J ( 5 0 '> ? ARLJC50) ?AT J(50) ?ARHJ<50) »ARQ(50> 
COMMON / B L K 2 / S K (33? 3 3 ) 

COMMON / BLK2 / S K ( 1 ? 1 ) 

COMMON / B L K 3 / S S T R ( 8 8 ) ? S S R 0 T (S3) 

COMMON 'BLKA/WR ( 8S ) ? Z ( 1 ? 1 ) ? PH I ( 88 38 ) ? OMEGA ( SS ) ? DYN ( 1 ? U ? 
C0MM0N/BLK6/WR ( 8 S )■.,-( i ;1 , put e 03 , qq ) , OMEGA ( 83 ; ? DYN ( L ? 1 ) ■> 
o SM '3 ? 1 > 

CDMM0M/BLK7/VLAM ( 3 > ? RL AM ( 3 ) ? CUEIB ( 20 ) 

C0MM0N/PLK3/L 

S' 0 M M ON / B L L 1 0 / S M i M V •' 1 ? 1 - 

COMMON XBLK 1 2/UREX 1 105 > • UIEX < 105 > 

0 0 M MON / RLE 13 'CMCQEI G ( 1 0 5 ) 

COMMON -'BLN 1 T/SENGD'. 1*1' 

COMMMN/BLM 4/SENSD ( 1 - 1 > 
r DMMOM ' T L ! ■ 3 A / A ’ 1?1) 

• ! 1 M M r ’ . ! / B L L 1 3 / P M P 1 

COMM 0 33 / B L K 2 0 / ] A U TO ? PEM A X 'AIM A X 

C0MM0N/BLK21/ST0M(33 ) 

COMMO N/BLL22/F ACT 

COMMON /BLK23/ISTATN 0(16) ? NST AT ■ CARR AY ( 1 6 ) 

C 0 M M 0 N / B L K 2 4 / R 1 V E C ( 1 6 ) 

C0MM0M/TI.L25/IME 
p n M M 0 I '/DLL 2 o / 1- H A T ( 10- 10' 

•_ ‘ , M M 0 N / B L L 2 7 / N TITLE •SO) 

.1 K I’ ) ¥ i’ * 

IMHilfUT DAI A SECT I IN 1 


f f T 'f r 'f 'f f " 7 f T I ? f 7 ? ' 

(±) f ' ! ‘ A il 3 0 ■ ( NT I 'X L ( I • - : - 1 - SO • 
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Appendix A- 4 


Program COPZ Documentation 

Damper Gain Optimization with Fixed Modal Damping Ratios 

I. INTRODUCTION 

COPZ is a damping optimization program to compute minimum 
damping gains for specified modal damping ratios using a diagonal 
damping matrix. Subroutine COMMIN is used in the optimization 
process . 

II. PROBLEM FORMULATIONS 

When a diagonal damping matrix [C] is introduced into an 
originally undamped system, the damped eigenvalue can be found by 
solving 

f(s) = det [ I + [fi(s)][6]] = 0 (A4.1) 

where 

s 

[S] 

[8] 


Let 

r t 

C 1 0 

• c 

0 N 

t c j 

The optimization problem is: 

find c. , i » 1 to N 
i c 
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= compelx eigenvalue of the damped system 
= submatrix of [AC] that contains only the non-zero 
terms of [AC] 

= corresponding sub-matrix of the receptance matrix R 


P P 
J « J t 

a n i j 

R = R = Z 

il J J «-l 2 2 

i j S + u> 

( 


N 

such that is minimized subject to the constraint 

f(S 4 ) =0, i = 1 to n e 


where 


S =-£(*> ±jVl”C 2 w 
l l n in 

i i 

To convert the above constrained optimization problem into one 
without constraints, a penalty function approach is used. The 
problem solved by COPZ can then be rephrased as 

find C . , i = 1 to N 
i c 

such that 


N N 

c E 

J - L c + Z W(£(S )) 
i = l i e = l i 


is minimized, where w is weighting function. Experience shows w = 
1000 yields good results. 

There axe N c + N £ unknowns in the optimized problems. The 
first N E of the unknowns are the unknown in (with 6 ^ 

specified). The last N c unknowns are the damper gains. 

II. INPUT TO COPZ 


A. Input from terminal (or input from TAPE 5). All are free format. 

1. Title. One line description of the problem. 

2. NDOF , NMODE , NE , NC 

NDOF = no. of system dof 
NMODE = no. of modes to be used 
NE = no. of specif icied £ 

NC = no. of dampers. 

3. JD( I ) , I = 1 , NC ) 

JD( I ) = Ith damper location 

4. X( I ) , I = 1 ,NC 
initial damper constant 

5 . NM( I ) , ETA( I ) , ( I = 1 to NE ) 

NM( I ) = mode no . 

ETA( I ) = C . = specified damping ratio 
WEIGHT 1 

WEIGHT = W is objective function J - 1000 usually 
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6. 


B . Modal Data Needed : ( same as routine SPSTGN ) 


Read from TAPE 12 : 

For I = 1 to NMODE 

WN( I ) , GM (Format 5X,2E16.9) 

(PHI( JJ,I),I=1,ND0F) (Format 5E16.9) 

Notes WN(I) = u)., GM( noticed), PHI(J,I) = fi. 
i 


IV- SAMPLE RUN OP COPZ 


Table 1 shows the job file to run COPZ, include data for a 
sample problem. For this sample problem, we place 6 dampers at dof 
1, 11, 39, 50, 78, and 88 of the 88-dof grillage model. We specify 
the model dampings to be 0.7 and 0.6 respectively for mode 4 and 5 
(the first vibration modes L The results are: 


C ± = 0.214, 

C 11 

= 0.268, 

C 39 

C 50 = °* 146 ' 

C 78 

= 0.207, 

n 

00 

00 


0.146 

0.267 


These results can be found in the end of Table 2. 


i 
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Table 1 CONTROL Cards for use of Program COPZ 
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Table 2 COPZ Sample Run Output 


i r ‘ O 

INPUT PROBLEM TITLE 

NASA GRILLAGE MODEL DAMPING OPTIMIZATION ? UAEY NMUDE 
I HPL'T : N DOF - HMGDE r ME .<• rJC 

o ~ n 

i M p :j p rt ^ m p r p LCC ^7 L n rl p 


INPUT INITIAL DAMPER CONST. uNT: 


INPUT PAIRS n - 


MODE MO, j DESIRED ZETA 


1 , ? 


NP'JT WEIGHT 


n m i - 1 1 


X H * ? X t % X T * * Y 

V 


t- * .1 i-, k\ I * | * * * * t * * t K 


I M C 0 N 




, nr; *. n c,> 


, l < r ' ! i ! *. r. • 


* 

■X 


i i 


< i 


X '■ : * * • t * 


s v i 1 
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Table 2 Continued 
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Table 2 Continued 
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Table 2 Continued 
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Table 2 Continued 
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Table 2 Continued 
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V. Listing of Program COPZ 
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Program Listing Continued 



Program Listing Continued 
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Program Listing Continued 
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Program Listing Continued 


r* 



* / T ’ 

1\ %. L ** L . 

• it-*/ i rnnjr / 2 ' £2 
r * pmjjt • ! £ T IC0UP7 
Ri! 


r i — t« 
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Appendix A- 5 


Program NORMAL Documentation 


Mode Shape Ortho Normalization 


I. INTRODUCTION 


NORMAL is a FORTRAN program that can be used to orthonormalize 
the mode shapes computed by EISPACK routines (such as RESV). For 
zero-frequency rigid body modes, Gram-Schmidt orthogonalization 
procedure ( see the Appendix A5 . 1 ) is used to make the rigid body 
modes orthogonal to each other with respect to the mass matrix. All 
modes are then normalized to unit generalized mass. Currently up to 
3 rigid body modes can be accepted. 

II. INPUT TO NORMAL 


The following data is needed to use NORMAL: 

1. Input N,NR,L where N = no. of dof of the system 

NR = no. of rigid body 
( <3 ) mode 

L = number of elastic modes to be normalized 

2. Read from TAPE 20 

For I = 1 to L 

WN(I), (V(J,I),J = 1,N) 

where 


WN( I ) = cj. = natural frequency of the ith mode, rad/sec 
V(J,I) = p.^ = mode shape coefficient at dof j of ith mode 


3 . Read from TAPE 10 


N J , ( AM( I ) , I = 1 , NJ ) 

NJ = no. of system dof 

where 

ith dof 

are created by unformatted WRITE in 


AM( I ) = m = lumped mass of 

Note: Both TAPE 10 and TAPE 20 

a prior analysis 
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III. PROGRAM LISTING 


F- F.OGF.AM :!GRM ( INF -JT . OUTPUT » TAF El 0 > TAFE12 * TAF £20 ) 


N-NO. OF THE SYSTEMS D . 0 . F . < 'JHNORMALIZED > » 

NF»- NO . OF RIGID EDDY MODES! 1-3) 

L -MO. OF ELASTIC MODES TO EE NORMALIZED TO UNIT GEN. MASS 
NATURAL FFEO. (RAD. /SEC. ) AND MODE SHAPES ARE STORED ON TAFi 
(USE r REE FORMAT) 

DIAGONAL ELEMENTS OF THE MASS MATRIX ARE STORED ON TAFEiO. 
THE NORMALIZED MODAL MATRIX IS WRITTEN ON TAPE! 2 . 


DIMENSION 'UN ( SS ( SR r ES '• » Z . EE f 3 ) » AM ( £?•> r PHI ( SS r BE ) » TT ( SS ) 


F' E A. I ‘f i : ! • r ! F, » L 
DO 200 1=1 >1 

PEAL" 20) W N •: I ) * ( '.*( Jt I ) » J = 1 iN ) 
read: io; j<j» •; AM<i)ri=irNj > 


NORMALIZE MODAL MATRIX 


DC 10C0 1 = 1 tt* 

G2 = 0.0 

DO 21 DO J-i»N 

G2 : 02 + V ( J r I : * V( 2 » I ) \ AM( J) 
6M = SORT-' 02) 

DO 2000 J = 1 • N 

V ( J • 1 ) -- V ( J » I ) / G M 

CONTINUE 


ORTHONCRMALIZE RIGID EDDY MODES 


IF (NR.LE.l) GO TO 9999 
DO 4100 I = 1 » N 
Z ( I » 1 ) = V ( I f 1 ) 

A = 0 

DO 4200 I=i?N 

A= A - V< I t i )* AM( I )tV( I»2) 

02 = 0 

DO 4200 1 = 1 rN 

2(1.2) = A * v ( I » 1 : r V ( I r 2 ) 

G2 = G2 -I- Z(If2)4ZCf2)*AM(I) 
GM = SORT (32) 

DO 4400 I = 1 r N 


•14 00 V ( I » 2 ) = Z ( I » 2 ) / 5 M 

IF (NR.LT.3) 0-0 TO ???? 


C2 

= 0 

1 

1 

Cl 

- 0 

t 

DO 

0100 

I = 1 > N 

W 

= Cl 

- Z(I»1 • tAM<I>*V<Xf3> 

C 2 

= C2 

- Z(If2)*AM(I)»V(I»3) 

G2 

= 0 

. 

I‘G 

C 2 00 

I = 1 » N i 

7 / T 

» 3 > = 

ci*Z(i.i) + c 2 1 z : r 2 ) i v ( i - z > 

u2 

= G2 

t Ztlr3)1Z(Ii3>*fsM(I) 

GM 

= SORT < G2 ) 

DC 

5300 

1 = 1 fN 


2 3 v 0 V < I * 3 > = Z ( I » 3 ) / G M 

PRINT * f WN(l)f (V(Jfl)t J«lrN) 
GM = 1.0 
1:0 6 000 I = 1 » N 
WRITE (12) UN ( I ) f GM 
cC-00 WRITE ( 12) ( V ( J » I ) > J = 1 » N ) 

COll FORMAT <5E 16. 9) 

6011 FORMAT (5Xr2E16. 9) 

? t 9 9 CONTINUE 

STOP 
END 
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APPENDIX A5.1 


Gram-Schmidt Orthogonal iz at ion Prodecure 

n independent vectors {v.} can be transformed into vectors 
{p. } orthogonal with respect £o a matrix [m] using the following 
recursion equations: 


where 


{p^ = (I ) 1 

{p) i = {r} i + « il {p} i + « i2 (P> 2 + 




(P) ^ [n>1 (r } • 
{ P ) T [m] {r }_. 
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APPENDIX B 


Active Damper Analysis - Linear 


Figure (B.l) shows system 2. 2. 2.1, it is essentially the same as 
Figure f5.6> except that definitions of displacement coordinates have 
been added. Note that the concept is essentially different to that 
used in the prototype design, which uses a long-stroke actuator 
without a linkage. 

Definitions 

C - Damping constant 

D - Voice coil, i.e., shaker, displacement 
f - Frequency 
F - shaker force 

K - Ratio of maximum structural amplitude to maximum voice-coil 

amplitude at design conditions 

m^ - shaker mass 

in - Movinq mass 

2 

P - Rate of energy absorption 

R - Linkage ratio 

W - Energy absorbed from structure 

x - Displacement of structure 

y - Displacement of m 2 

z - Relative displacement 

U) - Radial frequency 

subscripts 

AL - Amplitude limited 
CH - Characteristic 
C - Per cycle 
D - Design 
FL - Frequency limited 
M - Maximum 
0 - complex amplitude 

T - Total 

Referring to Figure (B.l), the basic damping equation is 

F/R = cdx/dt = m 2 d 2 y/dt 2 = m 2 (d 2 x/dt 2 + d 2 z/dt 2 ) 

and, assuming SHM with complex amplitudes of x Q , z q , respectively 

222 1/2 

Iz I = lx I{1 + (C /w m, )} ' 

O O 2 

If motion is amplitude-limited, so that 


I z 


o ' AL 


RD 


M 


then 


(B.l) 


(B . 2 ) 


(B . 3 ) 
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EQUIVALENT 
STRUCTURAL HASS 





whereas, if the motion is force-limited, we get, from B1 


' X o'fL = F M /wcR (B.5) 

Defining design conditions as those for which amplitude 
limits and force limits coincide, so that 


lx I nT = lx l__ 
o AL o FL 


'Vd 


and introducing the ratio 


(B.6) 


K ‘ Ix o' D /D M <B - 7 > 

one can eliminate c from ( B . 4 } and ( B . 5 } to get the biquadratic equation 

(R/K) 4 - (R/K) 2 - (f CH /Kf D ) 4 = 0 (B . 8 ) 

where f is the characteristic frequency given by 
CH 

f CH = 27 t * F M /m 2 D M } (B.9) 

and f, is the desiqn frequency. Solving A8 

<3 

(R/K) 2 = 1/2 + 1/2{1 + 4(f CH /Kf D ) 4 } 1/2 (B.10) 

thus-, assuming that K, f , and f are qiven, the linkaqe ratio R 

D CH, 

can be found, and the linkaqe can oe designed. The damping 
constant c can be found by substituting back into (B.4) and assuming 
desiqn conditions 

C = 277f D m 2 {(R 2 /K 2 ) - 1} 1/2 (B.ll) 

Note that the electronic system is designed so that the damping 
constant c is a constant regardless of frequency. Thus, on 
rearranging (B.4), and substituting from (B.ll) 


lx 

lx 


jsIfL 

I 

o 


RZK 

{(R/K 2 - l)(f 2 /f 2 ) + l } /2 


f < fr 


(B . 12) 


Equation (B.12) defines the ratio of structural amplitude at 
frequency f to the design amplitude. It is less then unity. For 
frequencies above the design frequency, the structural amplitude is 
force limited, so that, according to (B.l) 


f > f 


(B . 13 ) 


he 

lx 


I f 

O-FL _ _D . 

I f ' 

o D 


D 


The above two equations define structural amplitudes for 
which the damping system would remain linear under conditions of 
SHM. If the structure were to vibrate at a single frequency, they 
would constitute an envelope of permissible amplitudes. The 
situation with complicated structural responses would be much more 
difficult to analyze. 


The energy removed from the structure per cycle is 


W 


f 


cx dx = 27 T fclx. 


(B . 14) 


thus, at design conditions 

W D = 2rr 2 f D clx o l 2 D (B.15) 

and , on substitution from ( B . 5 ) and ( B.7 ) 


W = W (R/K) 
D m 


(B . 16 ) 


where 


W M " *Vk 

is the maximum energy which can be removed by the shaker. Since 
R/K is always greater than unity, W is always less than W^, the 
discrepancy being due to phase differences between the structural 
and damper motions. 

From ( B. 14 ) and (B.15) 

W _ 1x12 

_£ f Q 


CB.17) 


thus, for the amplitude limited case of (B.12) 


(f/f n ) (R 2 /K 2 ) 

D (RVK -1) (fp/f ) + 1 


while for the force-limited case of (B.13) 


W 

_C 



f > f 


D 


(B.19) 


(B.20) 
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The rate of removal of energy is given by 


P = f W 

c 

thus, above the design frequency 


= f D W D ; 


f > f T 


while, below the design frequency. 


P/P D 


2 2 2 2 
(f /f^) (R /K ) 

2 2 2 2 
(R /K -1) (fp/l ) 


we have from (B.19) 


+ 1 


Taking the Goodman V102 shaker as representative of current 
practice, we use the following values as typicals 

mass, m = 2 lb = 0,9072 kg 
1 m 

maximum force. F = 2 lb = 8.890 N 
M r 

P-P deflection, 2D = 0.2" = 0.00508 m 

H 

Max work per cycle, w (from (B.17)) = 0.0709 Joules 

N 

Max work per cycle per unit mass, W /m = 0.0782 J/kg 

H 


(B . 21) 


(B . 22) 


(B . 23 ) 


As an example, consider a design condition of one inch 
structural amplitude , I x_ I , at a frequency f of one Hertz . Then 
from (B.7 ) D 


- K - ,X o'd /D M = 10 

The following three sets of calculations are for three values of 
the movinq mass in 

2 

Case 12 3 Dimensions 


m 2 

0.09072 

0.9072 

9.072 

kg 

f ( from A9 ) 

31.26 

9.886 

3,126 

Hz 

R/K (from A10) 

3.207 

1.264 

1.005 


c ( from All ) 

1.737 

4.407 

5.707 

NS/m 

W D ( from A16 ) 

0.0221 

0.0561 

0.0705 

J 

“T * "l + ”2 

0.9979 

1.8144 

9.9792 

kg 

W D /m T 

0.0221 

0.0309 

0,0071 

J/kg 

P/m^, (from A22) 

0.0221 

0.0309 

0.0071 

W/kg 

It will be noted 

that case 2, for 

which the 

maonet and 

moving 


mass are equal, gives the best performance in terms of watts of 
energy removed per kg for frequencies above one Hz. The damping 
constant c achieved is also close to the maximum. If system 2. 2. 2.1 
in Figure (5.7) were used instead, the total mass m^ would be 0.0972 
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kg, and the performance of the system would otherwise be identical 
to case 2 above. Thus, values for w/m^ and P/m^ would be 0.0619 
J/kg and 0.619 w/kg respectively. This 1 value for W^/m^ is close to 
the maximum value of W^/m obtainable with this voice-coil shaker. 
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APPENDIX C 


Use of Piezoelectric Materials 


Polvvinylidene flouride (PVF 2 > exhibits the following 
relationship between strain e , normal electric field, E , and 
stress a (see Figure (C.l)) 1 3 


e. 


d 31 E 3 + 


(l/Cj 1 )o 1 


(C.l) 


where d 3l = transverse piezoeletric charge coefficient 
E 

C = modulus of elasticity 

Assuming SHM with maximum values of e , and a , such that 

e and a are 90° out of phase with each other, and 4^° out of 
phase with E 3 . 

°M - d 31 C ll E M /V * (C - 2) 

e M * d 31 E M /V2 (C - 3) 


and the maximum energy absorbed per unit volume is 

W M * 7T °H e M ' " d 31 C 11 E M 2/2 (C ‘ 4) 

while 

W M /m = W M /P = n d 31 C 11 E M 2/2p (C.5) 

Typical values are 

d 31 = 5 to 37 x 10~ 12 m/V 
= 1.6 to 3.8 x 10 9 N/m 2 
- 12 to 30 x 10 6 V/m 

p = 1.8 x 10 2 kg/m 2 
which gives the following range of values 

W /m = .005 to 4.1 J/Kg 

M 
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The most critical property is that of E^. Unless a 
considerable safety margin is allowed, there is a severe danger of 
dielectric breakdown in a space application, particular after a 
micrometeorite has damaged the material. 
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APPENDIX D 


Simulation of NASA Tests 


To help the actuator hardware design, a model of the prototype 
actuators will be put on the finite element model to simulater a 
vibration control test, initially, a linear simulation will be 
performed. The purpose is to determine the magnitude of the 
displacement and the actuator force under designed conditions. The 
assumptions of this simulation are 

( 1 ) the system is linear 

(2) the FEM of the NASA grillage (or beam) will be 
used as the undamped structure 

( 3 ) the actuator dynamics will be ignored . That is , 
the actuator will be assumed to be a perfect device 
capable of producing force f = cx under all conditions. 

( 4 ) the weight of each damper will be included as a 
discrete mass added to the system 

( 5 ) the excitation will be an initial displacement and/or 
initial velocity of the FEM. 

Formulations 


Let [M], [K] be the mass and stiffness of the undamped 
structure. Introducinq n dampers at dof , the damping 

matrix is c 1 

[c] - (Bn2][Bl T (D.l) 


where 


and 


[B] = [{e } I <e } 

1 2 


I )] 


nc 


nxn 



{e } = null vector except the jth element has a value of 1 

J 


In expanded form, 
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r 


0 


[C] = 

L 




J 


The contribution to the mass matrix due to the mass of the damper is 
[AM] = [B][&M][B] T 


or 


r 1 

l o I 


* Am 


[AM] = 

v 


* Am 


Am 


j 


where we assume all dampers are identical and have a mass of Am. 


The equation of motion of the system is 


[M]{x> + [ c ] { x } + [K] {x} = {0} 


where 


[M] = [M] + [AM] 
The initial conditions are 
{ x ( 0 ) } - { X Q } 

{i(0)} = {i Q } 


Equation (D.4) can be solved numerically using any of the numerical 
solution techniques, e.g. Wilson's 6 method or Newmark 0 method. 
Alternatively, a modal solution of (4) can be obtained by using the 
complex mode of the system (see Hurty and Rubinstein (Ref. 1 ) . 

Once the modes are known an analytical expression can be written 
for each mode in this approach. 

REFERENCE 


1. W.C. Hurty and M.F. Rubenstein, Dynamics of Structures , Prentice 
Inc. Englewood Cliffs, N.J. 1964 (ch. 9) . 
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(D. 2) 


(D . 3 ) 


(D.4) 


(D. 5 ) 


Kali, 


APPENDIX E 


1. Bounds on Transient Responses 
2. Time Domain Optimal Design 


1. Bounds on Transient Responses 


The disturbances acting on an Large Space Structure is mission 
dependent. However, one can say that these disturbances are usually 
of a transient type, e.g. firing of control thrusters for attitude 
control etc. The vibration control system must be designed to 
suppress these transient vibrations. One way of measuring the 
performance of a vibration control system is to evaluate how fast the 
control system can bring the response from some initial 
disturbed state to within an "acceptable” level. The response could 
be the displacement, velocity, or acceleration at certain 
"critical" locations of the LSS. For a large structure, it may not be 
easy to define which specific response is critical and the designer 
usually resorts to some aggregation of responses as the performance 
measure. One commonly used criterion is the minimization of 


V({Z}) = {Z} T [P]{Z}; z 


1 


X 

x 


j 


(E.l) 


where [P] is the weighting matrix. The estimation of transient 

behavior for the performance index can be derived according to 

Lyapunov's 2nd method, given a linear system governed by 

{z} = [A] {z} with initial state fz ) at t = t and where V({z},t ) 

o o o 

is given, the system attains a value of v given by 

V({z},t) =V({z },t )e ^ ^ in t-t seconds usinq Lyapunov’s 
o o o 

second method (Ref. 1,2). The following is a description of the 
use of Lyapunov’s second method in estimating transient responses 
(Ref. 6) and then indicate how it is applied to the control of LSS. 

Finally, some research issues will be raised. 

Given initial states (z } and the initial value V( {z },t ), 

o o o 

the system reaches the region 

V({z} f t) - V(x ,t )e~ V(t ~ t o ) (E.2) 

o o 


in t-t seconds, 
o 

For a linear time invariant system governed by the following 
state equations 
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{ z } = [A]{z> 


(E . 3 ) 


the parameter v in Eq. (E.2) is the minimum eigenvalue of the matrix 
tQl(P3 where [P] and [Q] are related by 

[A] T [P] + [PI CA] = -[Q] (E . 4 ) 

and [P] is the weighting matrix in the performance index, Eq. (E.l). 

Note that if [Q] is symmetric, so will be [p]. 

The above results can be applied to the control of LSS by 
letting 


[A] 


0 

1 

1 ! 

- 1 , 

i 

i 

i 

i i 

— < i 

i 

-M k 

-M c 

j 


(E . 5 ) 


where [A] is the dynamical matrix of the closed loop system. 

That is, in Eq. (E.5) the mass matrix [H] include the mass of the 
control devices; [kl is the stiffness matrix of the system plus the 
control law that is proportion to the displacement ( [x> } and [cl is 
the damping matrix which corresponds to the portion of control law 
that is proportional to rate ({x}), see Appendix F for detailed 
formulation. 


It should be noted that the use of Lyapunov’s second method 
requires that all the eigenvalues of the matrix [A] have negative 
real parts. For a free-free structure with no position feedback, 
the system matrix [A] has some zero eigenvalues. Thus, to use the 
above results directly we need to include positive feedback terms 
in the control law. 


The Lyapunov's second method provides a sound basis for 
transient performance estimates . To apply it to the design of 
vibration control of LSS, the following research needs to be 
carried out . 


(1) matrix [Q] must be positive definite, otherwise its 

eigenvalue could be zero or negative which renders the 
estimate (Eq. (E.2)) useless. This implies the performance 
index V({z),t) is restricted to quadratic form with 
weight matrix [P] that leads to a positive definite 
[Q]. Thus, [P] may not be specified arbitrarily. The 
question of how to pick matrix [P] needs further research. 
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(2) Alternatively, we can specify a positive definite [QJ and 
solve for [P] by Eq. (E.4>. This requires the solution of 
a set of n(n+l)/2 equation, where n is the order of matrix 
[A]. An efficient solution scheme of Eq. (E.41 needs to be 
develoepd ( or found in the literature ) . 

( 3 ) modify the formulation to allow it to accept LSS whose 
control law consists of rate feedback only. This may be 
achieved by formulating equations of motion for elastic 
modes only by filtering out the rigid body motion. 

( 4 ) application of the formulation in a truncated modal space . 
This will reduce the computational burden drastically. 

( 5 ) Formulate a synthesis procedure using the basic results . 

This may be achieved by iteractive methods or by a direct 
formulation. The results of the above researches can 
also be used in the next section. 

2. Time Domain Optimal Design 

Traditionally, time optimal design using optimal control 
theory with quadratic performance results in the solution of a higher 
order Riccatti equation, (Ref. 3). This will result in a control 
law with full state feedback which may be difficult to implement. 

For a vibration control system, it may be advisable to take another 
approach. That is one could fix the form of the control law then 
find the pairs that optimize some performance index. In general, 
nonlinear programming approaches can be used to find the optimal 
gains. This will be extremely costly since many transient 
responses are to be evaluated. If we limit ourself to a quadratic 
performance index, then the Lyapunov's second method can be used to 
evaluate the performance index when the system is driven by initial 
conditions. The general results (Ref. 6 ) will be stated first and 
then their application to LSS will be indicated. 

Consider the system 

{z} = [A] { z } 


where all eigenvalue sof [A] have negative real parts, or the 
original (0) is asympototically stable. The design problem is to 
adjust the elements in matrix [A] such that the performance index 

J - Jo fz} T [Q]fz)dt 

is minimized. It has been shown that 

J = J“ { z } T [ Q ] {z}dt - {z o ) T [P]{z o } 

where [Pi is the solution of the equation 


(E.6) 


(E . 7 1 
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(E.8) 


[A] T [P] + [P] [A] = -CQ] 

Note that [P] is symmetrical. The significance of Eq. (E.7) is that 
to evaluate the performance index, no transient analysis has to be 
performed at all. Rather, we just solve a set of n(n+l)/2 linear 
Eqs . (E.8) and then perform matrix multiplication . Since the 

elements of [P] is an implicit function of matrix [A] through Eq. 
(3), the following optimal procedure can be proposed: 

(1) from trial value of design parameters, form matrix [A] 

(2) solve [P] using Eq. (E.8) 

(3) evaluate J using Eq. (E.7) 

( 4 ) repeat the above process until J is minimized . 

When applied to LSS, we simply use the dynamic [A] as defined in 
Appendix F. The results of the research is used in the previous 
section is applicable here also. 


REFERENCES 

1. R.E. Kalman and J.E. Bertram, "Control System Analysis and 
Design Via the Second Method of Lyaponov" , Trans. ASME J. of 
Basic Engineering, 80, June 1960, pp. 371-400. 

2. K. Ogata, Modern Control Engineering , Prentice-Hall, Inc. 
Englewood Cliffs, New Jersey. 

3. Brian D. Anderson and John B. Moore, Linear Optimal Control , 
Prentice-Hall, Inc., Englewood Cliffs, New Jersey. 
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APPENDIX F 


Formulation of Dynamical Matrix for LSS with General Control Law 


Consider the actively controlled LSS 

[MJ[x] + [k Q ]{z} = [B ] { u } (F.l) 


where 


[k 1 is the stiffnes matrix of the uncontrolled structure (nxn) 

[M? is the mass matrix of the controlled system including the 
mass and inertia of the control devices ( sensors and 
actuators , etc ) ( nxn ) 

[B] is the control distribution matrix (NxN ) 

1 J u 

(u) is the control input (N xl) 

n 

Let the control law be 

{u} = - [ Ak ] { y } - [Ac]{y] (F.2) 


where 


{y} = [B d ]{x> 
{y} - [B u ] f x} 


displacement measurement (N,xl) 

d 

velocity measurement (N u xl) 


(F . 3 ) 


[AK] = displacement gain matrix (N u xN d ) 

[Acl = velocity qain matrix (N xN ) 

1 J • ■ u u 

[B^] = displacement measurement matrix (N^xN) 

[B u ] = velocity measurement matrix (N^xN) 

Substitute ( F . 3 ) into (F.2) 

(u) = -[Ak] [B d l fx] - [Ac][B u ]{i] (F.4) 

Substitute (F.4) into (F.l) 
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(F.5) 


[M3 { X } + [K o Hx} - -[BAkB d Hx} - [BAc BJfx) 
Equation (F.5) can be simplified to 

[M]{x> + [c] {x} + [k] {x} = {0} 


where 


[c] = [B][Ac][B v ] 

[K] = [K q ] + [AK] 

[AK] = [B][Ak][B d ] 

Note that in general, matrix [K] and [C] are arbitrary matrices, 
they may not be symmetric. 

Finally, we can put (F.6) in state space form as 


[z) = [A] {z} 


where 


[A] 


0 

1 

I 

-M _1 K 

-M~ c 

J 



1 


J 
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(F.6) 


(F.7) 


(F . 8 ) 


(F . 9 ) 


Appendix G 


Proposed Approaches for Actuator Failure 


Detection 

The actuator failure modes can be classified as: 

( a ) total failure : actuator does not produce a force in response 
to a command signal 

(b) random failure: actuator produces random force for a given 
input signal. 

The total failure mode may be considered in the design of 
electronics for the actuator. 

The random failure mode is more complex and may not be 
detected using hardware alone. In this case, a monitoring system 
to keep track of the performance of each actuator may be needed. 
Through the performance monitoring an adaptive control approach may 
be derived to detect andisolate failed actuators. Such an approach 
has been used in the detection of rate gyro failure in flight 
control systems [Ref. 1], The feasibility of this approach to our 
problem is being investigated. Literature search will be made to 
look for other approaches. 
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